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Abstract 

We introduce a near-linear complexity (geometric and meshless/algebraic) multi¬ 
grid/multiresolution method for PDEs with rough {L°°) coefficients with rigorous 
a-priori accuracy and performance estimates. The method is discovered through 
a decision/game theory formulation of the problems of (1) identifying restriction 
and interpolation operators (2) recovering a signal from incomplete measurements 
based on norm constraints on its image under a linear operator (3) gambling on 
the value of the solution of the PDE based on a hierarchy of nested measurements 
of its solution or source term. The resulting elementary gambles form a hierar¬ 
chy of (deterministic) basis functions of Hq{Q) (gamblets) that (1) are orthogonal 
across subscales/subbands with respect to the scalar product induced by the energy 
norm of the PDE (2) enable sparse compression of the solution space in i?o(D) (3) 
induce an orthogonal multiresolution operator decomposition. The operating dia¬ 
gram of the multigrid method is that of an inverted pyramid in which gamblets are 
computed locally (by virtue of their exponential decay), hierarchically (from fine 
to coarse scales) and the PDE is decomposed into a hierarchy of independent lin¬ 
ear systems with uniformly bounded condition numbers. The resulting algorithm is 
parallelizable both in space (via localization) and in bandwith/subscale (subscales 
can be computed independently from each other). Although the method is deter¬ 
ministic it has a natural Bayesian interpretation under the measure of probability 
emerging (as a mixed strategy) from the information game formulation and mul¬ 
tiresolution approximations form a martingale with respect to the hltration induced 
by the hierarchy of nested measurements. 


1 Introduction 

1.1 Scientific discovery as a decision theory problem 

The process of scientific discovery is oftentimes based on intuition, trial and error and 
plain guesswork. This paper is motivated by the question of the existence of a rational 
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decision framework that could be used to facilitate/guide this process, or turn it, to some 
degree, into an algorithm. In exploring this question, we will consider the problem of 
finding a method for solving (up to a pre-specified level of accuracy) PDEs with rough 
{L°°) coefficients as fast as possible with the following prototypical PDE (and its possible 
discretization over a fine mesh) as an example 

J — div (a(x)Vu(x)) = ^(x) x G g £ or g G 

= 0 on 

where fl is a bounded subset of (of arbitrary dimension d G N*) with piecewise 
Lipschitz boundary, a is a symmetric, uniformly elliptic d x d matrix with entries in 
L°°(fl) and such that for all x G and I G 

Amin(a)KP < l'^a{x)l < Amax(a)KP- (1-2) 

Although multigrid methods [41, 16, 50, 51, 101] are now well known as the fastest 
for solving elliptic boundary-problems and have successfully been generalized to other 
types of PDEs and computational problems [123], their convergence rate can be severely 
affected by the lack of regularity of the coefficients [37, 114]. Furthermore, although sig¬ 
nificant progress has been achieved in the development of multigrid methods that are, to 
some degree, robust with respect to meshsize and lack of smoothness (we refer in partic¬ 
ular to algebraic multigrid [94], multilevel finite element splitting [124], hierarchical basis 
multigrid [9, 24], multilevel preconditioning [106], stabilized hierarchical basis methods 
[107, 109, 110], energy minimization [65, 114, 122, 121, 108] and homogenization based 
methods [37, 34]), the design of multigrid methods that are provably robust with respect 
to rough {L°°) coefficients has remained an open problem of practical importance [17]. 

Alternative hierarchical strategies for the resolution of (1.1) are (1) wavelet based 
methods [18, 15, 3, 29, 38] (2) the Fast Multipole Method [49] and (3) Hierarchical ma¬ 
trices [52, 11]. Although methods based on (classical) wavelets achieve a multiresolution 
compression of the solution space of (1.1) in and although approximate wavelets 
and approximate L?' projections can stabilize hierarchical basis methods [109, 110], their 
applications to (1.1) are limited by the facts that (a) the underlying wavelets can per¬ 
form arbitrarily badly [7] in their 77 q(D) approximation of the solution space and (b) 
the operator (1.1) does not preserve the orthogonality between subscales/subbands with 
classical wavelets. The Fast Multipole Method and hierarchical matrices exploit the 
property that sub-matrices of the inverse discrete operator are low rank away from the 
diagonal. This low rank property can be rigorously proven for (1.1) (based on the ap¬ 
proximation of its Green’s function by sums of products of harmonic functions [10]) and 
leads to provable convergence (with rough coefficients), up to the pre-specified level of 
accuracy e in L^-norm, in 0{N\n^ Aj operations ([10] and [11, Thm. 2.33 and 

Thm. 4.28]). Can the problem of finding a fast solver for (1.1) be, to some degree, 
reformulated as an Uncertainty Quantification/Decision Theory problem that could, to 
some degree, be solved as such in an automated fashion? Can discovery be computed? 
Although these questions may seem unorthodox their answer appears to be positive: this 
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paper shows that this reformulation is possible and leads to a multigrid/multiresolution 
method/algorithm solving (1.1), up to the pre-specihed level of accuracy e in ff^-norm 
(i.e. finding u^'PP such that Hii — tt’^PP||j:^i(Q) < e||9||H-i(n) for an arbitrary g decomposed 
over N degrees of freedom), in C)(iV In^'^ (max(^, operations (for e ~ 

the hierarchical matrix method achieves e-accuracy in norm in 0{N N) op¬ 

erations and the proposed multiresolution method achieves e-accuracy in norm in 
0(A^In^'^ A^) operations). For subsequent solves (i.e. if (1.1) needs to be solved for 
more than one g) then the proposed multiresolution method achieves accuracy e ~ N~d 
in Ff^-norm in 0(A^ln'^'''^ A^) operations (we refer to Subsection 5.4 and in particular 
to Table 1 for a detailed complexity analysis of the proposed method, which can also 
achieve sublinear complexity if one only requires L^-approximations). 

The core mechanism supporting the complexity of the method presented here is the 
fast decomposition of into a direct sum of linear subspaces that are orthogonal 

(or near-orthogonal) with respect to the energy scalar product and over which (1.1) 
has uniformly bounded condition numbers. It is, to some degree, surprising that this 
decomposition can be achieved in near linear complexity and not in the complexity of 
an eigenspace decomposition. Naturally [86], this decomposition can be applied to the 
fast simulation of the wave and parabolic equations associated to (1.1) or to its fast 
diagonalization. 

The essential step behind the automation of the discovery/design of scalable numer¬ 
ical solvers is the observation that fast computation requires repeated computation with 
partial information (and limited resources) over hierarchies of levels of complexity and 
the reformulation of this process as that of playing underlying hierarchies of adversarial 
information games [111, 112]. 

Although the problem of finding a fast solver for (1.1) may appear disconnected 
from that of finding statistical estimators or making decisions from data sampled from an 
underlying unknown probability distribution, the proposed game theoretic reformulation 
is, to some degree, analogous to the one developed in Wald’s Decision Theory [113], 
evidently influenced by Von Neumann’s Game Theory [111, 112] (the generalization of 
worst case Uncertainty Quantihcation analysis [83] to sample data/model uncertainty 
requires an analogous game theoretic formulation [80], see also [79] for how the underlying 
calculus could be used to guide the discovery of new Selberg identities). We also refer 
to subsection 1.3 for a review of the correspondence between statistical inference and 
numerical approximation. 

1.2 Outline of the paper 

The essential difficulty in generalizing the multigrid concept to PDFs with rough coeffi¬ 
cients lies in the fact that the interpolation (downscaling) and restriction (upscaling) op¬ 
erators are, a priori, unknown. Indeed, in this situation, piecewise linear finite-elements 
can perform arbitrarily badly [7] and the design of the interpolation operator requires 
the identification of accurate basis elements adapted to the microstructure a(x). 

This identification problem has also been the essential difficulty in numerical homoge- 
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nization [117, 6 , 4, 18, 59, 33, 84, 17]. Although inspired by classical homogenization ideas 
and concepts (such as oscillating test functions [ 68 , 36, 35], cell problems/correctors and 
effective coefficients [13, 90, 1, 72, 39, 46], harmonic coordinates [61, 6 , 4, 76, 12, 2, 84], 
compactness by compensation [100, 45, 67, 14]) an essential goal of numerical homoge¬ 
nization has been the numerical approximation of the solution space of ( 1 . 1 ) with arbi¬ 
trary rough coefficients [84], i.e., in particular, without the assumptions found in classical 
homogenization, such as scale separation, ergodicity at fine scales and e-sequences of op¬ 
erators (otherwise the resulting method could lack robustness to rough coefficients, even 
under the assumption that coefficients are stationary [ 8 ]). Furthermore, to envisage ap¬ 
plications to multigrid methods, the computation of these basis functions must also be 
provably localized [5, 85, 64, 48, 87, 58] and compatible with nesting strategies [87]. In 
[77] , it has been shown that this process of identihcation (of accurate basis elements for 
numerical homogenization), could, in principle, be guided through its reformulation as a 
Bayesian Inference problem in which the source term (7 in (1.1) is replaced by noise ^ and 
one tries to estimate the value of the solution at a given point based on a finite number 
of observations. In particular it was found that Rough Polyharmonic Splines [87] and 
Polyharmonic Splines [54, 30, 31, 32] can be re-discovered as solutions of Gaussian filter¬ 
ing problems. This paper is inspired by the suggestion that this link between numerical 
homogenization and Bayesian Inference (and the link between Numerical Quadrature 
and Bayesian Inference [92, 27, 97, 74, 75]) are not coincidences but particular instances 
of mixed strategies for underlying information games and that optimal or near optimal 
methods could be obtained by identifying such games and their optimal strategies. 

The process of identification of these games starts with the (Information Based Com¬ 
plexity [119]) notion that computation can only be done with partial information. For 
instance, since the operator ( 1 . 1 ) is infinite dimensional, one cannot directly compute 
with u G Hq (II) but only with finite-dimensional features of u. An example of such finite¬ 
dimensional features is the m-dimensional vector Um ■= (Jq ■ ■ ■, Jq ufim) obtained 
by integrating the solution u of ( 1 . 1 ) against m test/measurement functions (pi G L^(II). 
However to achieve an accurate approximation of u through computation with Um one 
must fill the information gap between Um and u (i.e. construct an interpolation operator 
giving u as a function of Um)- We will, therefore, reformulate the identification of this 
interpolation operator as a non-cooperative (min max) game where Player I chooses the 
source term g (1.1) in an admissible set/class (e.g. the unit ball of L^(H)) and Player II 
is shown Um and must approximate u from these incomplete measurements. Using the 
energy norm 

ll^lla {x)a{x)Vu{x) dx, (1-3) 

Jn 

to quantify the accuracy of the recovery and calling u* Player I’s bet (on the value of 
u), the objective of Player I is to maximize the approximation error ||tt — M*||a, while 
the objective of Player II is to minimize it. A remarkable result from Game Theory 
(as developed by Von Neumann [111], Von Neumann and Morgenstern [112] and Nash 
[70]) is that optimal strategies for deterministic zero sum finite games are mixed (i.e. 
randomized) strategies. Although the information game described above is zero sum, 
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it is not finite. Nevertheless, as in Wald’s Decision Theory [113], under sufficient regu¬ 
larity conditions it can be made compact and therefore approximable by a finite game. 
Therefore although the information game described above is purely deterministic (and 
has no a priori connection to statistical estimation), under compactness (and continuity 
of the loss function), the best strategy for Player I is to play at random by placing a 
probability distribution tt/ on the set of candidates for g (and select g' as a sample from 
TT/) and the optimal strategy for Player II is to place a probability distribution vr// on 
the set of candidates for g and approximate the solution of (1.1) by the expectation of 
u (under tt// used as a prior distribution) conditioned on the measurements ucjji. 

Although the estimator employed by Player II may be called Bayesian, the game 
described here is not (i.e. the choice of Player I might be distinct from that of Player II) 
and Player II must solve a min max optimization problem over ttj and ttjj to identify an 
optimal prior distribution for the Bayesian estimator (a careful choice of the prior also 
appears to be important due to the possible high sensitivity of posterior distributions 
[81, 79, 82]). Although solving the min max problem over ttj and ttjj may be one way of 
determining the strategy of Player II, it will not be the method employed here. We will 
instead analyze the error of Player IPs approximation as a function of Player IPs prior 
and the source term g picked by Player I. Furthermore, to preserve the linearity of the 
calculations we will restrict Player IPs decision space (the set of possible priors vr//) to 
Gaussian priors on the source term g. Since the resulting analysis is independent of the 
structure of (1.1) and solely depends on its linearity we will first perform this investi¬ 
gation, in Section 2, in the algebraic framework of linear systems of equations, identify 
Player IPs optimal mixed strategy and show that it is characterized by deterministic 
optimal recovery and accuracy properties. The mixed strategy identified in Section 2 
will then be applied in 3 to the numerical homogenization of (1.1) and the discovery of 
interpolation interpolators. In particular, it will be shown that the resulting elementary 
gambles form a set of deterministic basis functions (gamblets) characterized by (1) op¬ 
timal recovery and accuracy properties (2) exponential decay (enabling their localized 
computation) (3) robustness to high contrast. 

To compute fast, the game presented above must not be limited to filling the infor¬ 
mation gap between Um £ and u G PIq(D). This game must be played (and repeated) 
over hierarchies of levels of complexity (e.g. one must fill information gaps between 
and then and etc...). We will therefore, in Section 4, consider the (hier¬ 
archical) game where Player I chooses the r.h.s of (1.1) and Player II must (iteratively) 
gamble on the value of its solution based on a hierarchy of nested measurements of u 
(from coarse to fine measurements). Under Player IPs mixed strategy (identified in Sec¬ 
tion 2 and used in Section 3), the resulting sequence of multi-resolution approximations 
forms a martingale. Conditioning and the independence of martingale increments lead 
to the hierarchy of nested interpolation operators and to the multiresolution orthogonal 
decomposition of (1.1) into independent linear systems of uniformly bounded condition 
numbers. The resulting elementary gambles (gamblets) (1) form a hierarchy of nested 
basis functions leading to the orthogonal decomposition (in the scalar product of the 
energy norm) of P7 q(D) (2 ) enable the sparse compression of the solution space of (1.1) 
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(3) can be computed and stored in near-linear complexity by solving a nesting of linear 
systems with uniformly bounded condition numbers (4) enable the computation of the 
solution of ( 1 . 1 ) (or its hyperbolic or parabolic analogues) in near-linear complexity. 
The implementation and complexity of the algorithm are discussed in Section 5 with 
numerical illustrations. 

1.3 On the correspondence between statistical inference and numerical 
approximation 

As exposed by Diaconis [27], the investigation of the correspondence between statisti¬ 
cal inference and numerical approximation can be traced back to Poincare’s course in 
Probability Theory [92]. It is useful to recall Diaconis’ compelling example [27] as an 
illustration of this conection. Let / : [0,1] —)• M be a given function and assume that we 
are interested in the numerical approximation of fg f{t) dt. The Bayesian approach to 
this quadrature problem is to (1) Put a prior (probability distribution) on continuous 
functions C[0,1] (2) Calculate / at xi, X 2 , ■ ■ ■ ,Xn (to obtain the data (/(xi),..., f{xn))) 
(3) Compute a posterior (4) Estimate fg f{t) dt by the Bayes rule. If the prior on C[0,1] 
is that of a Brownian Motion (i.e. f{t) = Bt where Bt is a Brownian motion and Bq is 
normal), then E[/(x)|/(xi), ..., f{xn)] is the piecewise linear interpolation of / between 
the points xi,... ,x„ and one re-discovers the trapezoidal quadrature rule. If the prior 
on C[0,1] is that of the first integral of a Brownian Motion (i.e. f{t) fg Bg ds) then 
the posterior E[/(x)|/(xi), .. ., /(x„)] is the cubic spline interpolant and integrating k 
times yields splines of order 2A: -|- I. 

Subsequent to Poincare’s early discovery [92], Sul’din [102] and (in particular) Larkin 
[62] initiated the systematic investigation of the correspondence between conditioning 
Gaussian measures/processes and numerical approximation. As noted by Larkin [62], de¬ 
spite Sard’s introduction of probabilistic concepts in the theory of linear approximation 
[95] , and Kimeldorf and Wahba’s exposition [60] of the correspondence between Bayesian 
estimation and spline smoothing/interpolation, the application of probabilistic concepts 
and techniques to numerical integration/approximation “attracted little attention among 
numerical analysts” (perhaps due to the counterintuitive nature of the process of random¬ 
izing a known function). However, a natural framework for understanding this process 
of randomization can be found in the pioneering works of Wozniakowski [118], Packel 
[ 88 ], and Traub, Wasilkowski, and Wozniakowski [104] on Information Based Complex¬ 
ity [71, 119], the branch of computational complexity that studies the complexity of 
approximating continuous mathematical operations with discrete and finite ones up a to 
specified level of accuracy. Indeed the concept that numerical implementation requires 
computation with partial information and limited resources emerges naturally from In¬ 
formation Based Complexity, where it is also augmented by concepts of contaminated 
and priced information associated with, for example, truncation errors and the cost of 
numerical operations. In this framework, the performance of an algorithm operating on 
incomplete information can be analysed in the usual worst case setting or the average 
case (randomized) setting [93, 73] with respect to the missing information. Although 
the measure of probability (on the solution space) employed in the average case setting 
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may be arbitrary, as observed by Packet [88], if that measure is chosen carefully (as the 
solution of a game theoretic problem) then the average case setting can be interpreted 
as lifting a (worst case) min max problem (where saddle points of pure strategies do 
not, in general, exist) to a min max problem over mixed (randomized) strategies (where 
saddle points do exist [111, 112]). As exposed by Diaconis [27] (see also Shaw [97]) the 
randomized setting also establishes a correspondence between Numerical Analysis and 
Bayesian Inference providing a natural framework for the statistical description of nu¬ 
merical errors (in which confidence intervals can be derived from posterior distributions). 
Furthermore [89, 27], classical min max numerical quadrature rules can be formulated 
as solutions of Bayesian inference problems with carefully chosen priors [27] and, as 
shown by O’Hagan [74, 75], this correspondence can be exploited to discover new and 
useful numerical quadratures rules. As envisioned by Skilling [99], by placing a (care¬ 
fully chosen) probability distribution on the solution space of an ODE and conditioning 
on quadrature points, one obtains a posterior distribution on the solution whose mean 
may coincide with classical numerical integrators such as Runge-Kutta methods [96] . As 
shown in [23] the statistical approach is particularly well suited for chaotic dynamical 
systems for which deterministic worst case error bounds may provide little information. 
While in [99, 96, 23] the probability distribution is directly placed on the solutions space, 
for PDEs [77] argues that the prior distribution must be placed on source terms (or on 
the image space of an integro-differential operator) and propagated/filtered through the 
inverse operator to reflect the structure of the solution space. In particular [77] shows 
that this process of filtering noise with the inverse operator, when combined with con¬ 
ditioning, produces accurate finite-element basis functions for the solution space whose 
deterministic worst case errors can be bounded by standard deviation errors using the 
reproducing kernel structure of the covariance function of the filtered Gaussian field. As 
already witnessed in [23, 96, 77, 56, 55, 19, 25], it is natural to expect that the possibili¬ 
ties offered by combining numerical uncertainties/errors with model uncertainties/errors 
in a unified framework will stimulate a resurgence of the statistical inference approach 
to numerical analysis. 

2 Linear Algebra with incomplete information 

2.1 The recovery problem 

The problem of identifying interpolation operators for (1.1) is equivalent (after discretiza¬ 
tion or in the algebraic setting) to that of recovering or approximating the solution of 
a linear system of equations from an incomplete set of measurements (coarse variables) 
given known norm constraints on the image of the solution. 

Let n > 2 and A be a known real invertible n x n matrix. Let b be an unknown 
element of M"". Our purpose is to approximate the solution x of 

Ax = b (2.1) 
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based on the information that ( 1 ) a: solves 

= y, (2.2) 

where (the measurement matrix) is a known, rank m, m x n real matrix such that 
m < n and y (the measurement vector) is a known vector of M™, and ( 2 ) the norm 
h'^T-^h of b is known or bounded by a known constant (e.g., ^b < 1), where T ^ 
is a known positive definite n x n matrix (with T~^ being the identity matrix as a 
prototypical example). Observe that since m < n, the measurements ( 2 . 2 ) are, a priori, 
not sufficient to recover the exact value x. 

As described in Section 1 , by formulating this recovery problem as a (non-cooperative) 
information game (where Player I chooses b and Player II chooses an approximation x* 
of X based on the observation <hx), one (Player II) is naturally lead to search for mixed 
strategy in the Bayesian class by placing a prior distribution on b. The purpose of this 
section is to analyze the resulting approximation error and select the prior distribution 
accordingly. To preserve the linearity (i.e. simplicity and computational efficiency) of 
calculations we will restrict Player IPs decision space to Gaussian priors. 

2.2 Player I’s mixed strategy 

We will therefore, in the first step of the analysis, replace b in (2.1) by a centered 
Gaussian vector of M"" with covariance matrix Q (which may be distinct from T) and 
consider the following stochastic linear system 

AX = ^. (2.3) 

The Bayesian answer (a mixed strategy for Player II) to the recovery problem of Section 
2 is to approximate x by the conditional expectation E[W|<I>A = y]. 

Theorem 2.1. The solution X of (2.3) is a centered Gaussian vector o/M'^ with co- 
variance matrix 

K = A-^Q{A-^f. (2.4) 

Furthermore, X conditioned on the value ^X = y is a Gaussian vector ofMA with mean 
E[X|4'A = y] = Ty, and of covariance matrix , where T is the n x m matrix 

^ ■= (2.5) 

and is the rank n — m positive nxn symmetric matrix defined by := K — 

Proof. (2.4) simply follows from X = Since X is a Gaussian vector, E[X|4>X = 

y] = 'ky where T is a n x m matrix minimizing the mean squared error E [|X — Md^Xp] 
over all nxm matrices M. We have E[|X —M4>Xp] = Trace[X]-|-Trace[M‘kX<I>^M^] — 
2 Trace[‘I>XM] whose minimum is achieved for M = T as defined by (2.5). The co- 
variance matrix of X given ‘kX = y is then obtained by observing that for v G M”, 
v'^K'^v = E [\v'^X - u'^T 4 >Xp] = v'^Kv - u^T4>Xu. □ 
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2.3 Variational/optimal recovery properties and approximation error 

For a n X n symmetric positive definite matrix M let (t)^ be the (scalar) product 

on M"" defined by: for u,v G M”, (u, := u^Mv and write ||f||M := 

corresponding norm. When M is the identity matrix then we write (u, n) and HuH the 
corresponding scalar product and norm. For a linear subspace V of MT' we write Pv,m 
for the orthogonal projections onto V with respect to the scalar product For a 

(possibly rectangular) matrix B we write Im(i?) the image (range) of B and Ker(S) the 
null space of B. For an integer n, let /„ be the n x n identity matrix. 

Theorem 2.2. For w G M™, is the unique minimizer of the following quadratic 
problem 

{ Minimize , 

\ ’ //f-i (2.6) 

Subject to = w and u G M" . 

In particular, v = '^y, the Bayesian approximation of the solution of (2.1), is the unique 
minimizer of ||j4u||q-i under the measurement constraints ‘hu = y. Furthermore, it also 
holds true that (1) = Im (2) Im('I^) is the orthogonal complement o/Ker(‘h) with 

respect to the product (•, •)^-i and (3) 'I'<1> = and In — = Fker($),ii'-i • 

Proof. First observe that (2.5) implies that ‘h'l/ = I^ where Im is the identity mxm ma¬ 
trix. Therefore $('I't(;) = w. Note that (2.5) implies that for all 2 ; G M"*, = 

2 ;^(<l>iir<h^) ^<l>u. Therefore if u G Ker(^>) then = 0 for all 2 ; G M”^. Con¬ 
versely if = 0 for all 2 ; G M”* then v must belong to Ker(<h). Since the 

dimension of Im('I') is m and that of Ker(<l>) is n — m we conclude that Im(^') is the 
orthogonal complement Ker(<h) with respect to the product (•, •)^_i and in particular, 

('I'rc, = 0, yw G and Vu G M” such that = 0 . (2.7) 

Let w G and u G M” such that = w. Since 'Lrc — u G Ker(<l>), it follows from (2.7) 
that + (^v — 'I>w,v — Therefore 'Lu; is the unique 

minimizer of over all v G MF such that ‘hu = w. Now consider / G M”, since 

Im('I') = Im(iL<L^) and Im(^') is the orthogonal complement of Ker(<l>) with respect to 
the product there exists a unique 2 ; G and a unique g G Ker(^>) such that 

/ = z -|- g. Since 'I'<h = iL<l>^(<l>7L<h^)“^<l>, it follows that = K^^z and (/„ — 
'^^)f = g. We conclude by observing that g = PKer{^),K-^ /• □ 

Theorem 2.3. For v G M”', w* = d>u is the unique minimizer of ||u — over 

all w G M”^. In particular, ||u — = min^gigm ||ti — iL<L^ 2 ;||^-i and if x is the 

solution of the original equation (2.1), then ||x — "I>y\\x-^ = niin^gRm ||x — = 

minjjgRm ||x — iL<l>^2;||;^-i. 

Proof. The proof follows by observing that v — belongs to the null space of $ which, 
from Theorem 2.2, is the orthogonal complement of the image of 'L with respect to the 
scalar product defining the norm || • Observe also that the image of ^ is equal to 

that of □ 
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Remark 2.4. Observe that, from Theorem 2.2, v — 4^$?; spans the null space o/$, and 
= II?; — 4'$?;||^_i + ||4'4>?;||^_i. Therefore if D is a symmetric positive definite 
n X n matrix then sup^gj^n ||?; — 'I'4'?;||^/||?;||^-i = sup^gj|n^,j,^=o Ikll^’/ll^llE'-i • P®’"" 

ticular, if x is the solution of ( 2 . 1 ) and y the vector in ( 2 . 2 ), then ||x — 4'y||^/||6||Q-i < 
sup^g]Rn_$.u=o Ikllo/ll^lli^-i the right hand side is the smallest constant for which 
the inequality holds (for all b). 


Remark 2.5. A simple calculation (based on the reproducing Kernel property (^v, K.^f^ 

Vi) shows that if x is the solution of ( 2 . 1 ) and y the vector in ( 2 . 2 ), then 

{x — ^y)i < {Kf(f^\\b\\Q-i, i.e. the variance of the ith entry of the solution of the 
stochastic system (2.3) conditioned on = y, controls the accuracy of the approxima¬ 
tion of the ith entry of the solution of the deterministic system (2.1). In that sense, the 
role of is analogous to that of the power function in radial basis function interpolation 
[116, 40 ] and that of the Kriging function [120] in geostatistics [69]. 


2.4 Energy norm estimates and selection of the prior 

We will from now on assume that A is symmetric positive definite. Observe that in this 
situation the energy norm || • ^ is of practical significance for quantifying the approxi¬ 
mation error and Theorem 2.3 leads to the estimate \\x — ^y\\j^-i = min^gRm \\Q~ 2 b — 
Q~ 2 A 2 K 2 ^'^z\\ which simplifies to the energy norm estimate expressed by Corollary 
(2.6) under the choice Q = A (note that K~^ = A under that choice). 

Corollary 2.6. If A is symmetric positive definite and Q = A, then for v G 'MT, 
llu — = minzgRm Hu — z\\a- Therefore, if x is the solution of (2.1) and y 

the vector in (2.2), then \\x — 'kym = min^gRm ||x — 'krum = min^gKm \\x — z\\a- 

In particular 

llx — = min \\A~2b — A~2^'^z\\. ( 2 . 8 ) 

Remark 2.7. Therefore, according to Corollary 2.6, ifQ = A, then Ty is the Galerkin 
approximation of x, i.e. the best approximation of x in || • \\A-norm in the image of 
'k (which is equal to the image of ). This is interesting because ^'y is obtained 

without the prior knowledge of b. 

Corollary 2.6 and Remark 2.7 motivate us to select Q = A as the covariance matrix 
of the Gaussian prior distribution (mixed strategy of Player II). 


2.5 Impact and selection of the measnrement matrix $ 

It is natural to wonder how good this recovery strategy is (under the choice Q = A) 
compared to the best possible function of y and how the approximation error is impacted 
by the measurement matrix <k. If the energy norm is used to quantify accuracy, then 
the recovery problem can be expressed as finding the function 0 of the measurements 
y minimizing the (worst case) approximation error infg) sup||ft||<]^ \\x — 6l(y)m/||&|| with 
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X = A~^h and y = ^A~^b. Writing 0 < Ai(A) < • • • < An(A), the eigenvalues of A 
in increasing order, and ai,..., On, the corresponding eigenvectors, it is easy to obtain 
that (1) the best choice for would correspond to measuring the projection of x on 
spanjai,..., am} and would lead to the worst approximation error 1 /Am+i and (2) the 
worst choice would correspond to measuring the projection of x on a subspace orthogonal 
to a\ and would lead to the worst approximation error l/\/Ai. Under the decision Q = A 
the minimal value of (2.8) is also 1/y^Am+i and achieved for Im(<h^) = spanjoi,..., am} 
and the maximal value of (2.8) is l/\/Ai and achieved when Im(<h^) is orthogonal to ai. 
The following theorem, which is a direct application of (2.8) and the estimate derived in 
[53, p. 10] (see also [66]), shows that, the subset of measurement matrices that are not 
nearly optimal is of small measure if the rows of are sampled independently on the 
unit sphere of M”. 

Theorem 2.8. If ^ is a n x m matrix with i.i.d. AA(0,1) (Gaussian) entries, Q = A, 
X is the solution of the original equation (2.1), and 2 < p then with probability at least 
1 - 3p~P, ||x - Ty||A/||?)|| < (1 + 9y/m+py/n)/^yXm+l■ 

Although the randomization of the measurement matrix [42, 63, 43, 78] can be an 
efficient strategy in compressed sensing [105, 21, 20, 28, 44, 22] and in Singular Value 
Decomposition/Low Rank approximation [53], we will not use this strategy here because 
the design of the interpolation operator presents the (added) difficulty of approximating 
the eigenvectors associated with the smallest eigenvalues of A rather than those asso¬ 
ciated with the largest ones. Furthermore, T has to be computed efficiently and the 
dependence of the approximation constant in Theorem 2.8 on n and m can be prob¬ 
lematic if sharp convergence estimates are to be obtained. We will instead select the 
measurement matrix based on the transfer property introduced in [14] and given in a 
discrete context in the following theorem. 

Theorem 2.9. If A is symmetric positive definite, Q = A and x is the solution of the 
original equation (2.1), then for any symmetric positive definite matrix B, we have 


inf 


v'^Bv 


v'^Av 5:eK™ 


min ||6 — 


B- 


< ||x —Tym < sup 


v'^Bv 


v'^Av srSM'" 


min \\b — ^'^z\\b-i 


(2.9) 


Proof. Corollary 2.6 implies that if x is the solution of the original equation (2.1), then 
\\x — Tym = min^gffim ||5 — 4>'^2;||^-i. We finish the proof by observing that if A and B 
are symmetric positive definite matrices such that aiB < A < a 2 B for some constants 
01,02 > 0 then < A~^ < af^B~^. □ 

Therefore according to Theorem 2.9, once a good measurement matrix <1> has been 
identified for a symmetric positive definite matrix B such that aiB < A, the same 
measurement matrix can be used for A at the cost of an increase of the bound on the 
error by the multiplicative factor o^ ' . As a prototypical example, one may consider a 
(stiffness) matrix A obtained from a finite element discretization of the PDF (1.1) and 
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B may be the stiffness matrix of the finite element discretization of the Laplace Dirichlet 
PDE 

— Au'{x) = g{x) on Q with u' = 0 on dO., (2-10) 

obtained from the same finite-elements (e.g. piecewise-linear nodal basis functions over 
the same fine mesh Th)- Using the energy norm (1.3), Theorem 2.9 and Remark 2.7 
imply the following proposition 

Proposition 2.10. Letuh (resp. he the finite element approximation of the solution 
u of (1.1) (resp. the solution u' of (2.10)J over the finite nodal elements ofTh- Let uh 
( resp. u'^) he the finite element approximation of the solution u of (1.1) (resp. the 
solution u' of (2.10)j over linear space spanned hy the rows of (resp. over the 

linear space spanned hy the rows of ). It holds true that 

—===\\Uh-Ujj\\H^(n) < \\uh-UH\\a < / . ~ (2-11) 

Observe that the right hand side of (2.11) does not depend on A ma x(a). therefore if 
Amin (a) = 1, then the error bound on \\uh — unWa does not depend on the contrast of a 
(i.6. Amax (a)/Amin (o.)). 


3 Numerical homogenization and design of the interpola¬ 
tion operator in the continuous case 

We will now generalize the results and continue the analysis of Section 2 in the con¬ 
tinuous case and design the interpolation operator for (1.1) in the context of numerical 
homogenization. 

3.1 Information Game and Gamblets 

As in Section 2 we will identify the interpolation operator (that will be used for the 
multigrid algorithm) through a non cooperative game formulation where Player I chooses 
the source term g (1.1) and Player II tries to approximate the solution u of (l.I) based 
on a finite number of measurements (/^ U()>i)i<i<m obtained from linearly independent 
test functions (fi G L^(U). As in Section 2, this game formulation, motivates the search 
for a mixed strategy for Player II that can be expressed by replacing the source term g 
with noise We will therefore consider the following SPDE 

f-div (^a(x)Vu(x)) = ^(x) x G U; 

= 0 on dll, 

where U and a are the domain and conductivity of (1.1). As in Section 2, to preserve the 
computational efficiency of the interpolation operator we will assume that ^ is a centered 
Gaussian field on U. The decision space of Player II is therefore the covariance function 
of Write C the differential operator — div(aV) with zero Dirichlet boundary condition 
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mapping Hq{Q) onto Motivated by the analysis (Remark 2.7) of Subsection 2.4 

(which can be reproduced in the continuous case) we will select the covariance function 
of ^ (Player IPs decision) to be C. Therefore, under that choice, for all / G 

f(x)^(x) dx is a Gaussian random variable with mean 0 and variance fCf = ||/||q 
where ||/||a is the energy norm of / defined in (1.3). Introducing the scalar product on 
defined by 

{v,w) := [ {VvfaVw, (3.2) 

Jn 

recall that if (ei, 62 ,...) is an orthonormal basis of (77 q (fl), || • ||a) diagonalizing C, then ^ 
can formally be represented as ^ (where the Xi are i.i.d. A/’(0,1) random 

variables) and, therefore, ^ can also be identified as the linear isometry mapping Hq{Q) 
onto a Gaussian space and / = {f, ei)^ei onto f{x)({x) dx = (/, 

Observe also that [77], if is White Noise on 0 (i.e. a Gaussian field with covariance 
function 6{x — y)) then ^ can be represented as ^ . Furthermore [77, Prop. 3.1] 

the solution of (3.1) is Gaussian field with covariance function G{x,y) (where G is the 
Green’s function of the PDF (1.1), i.e. CG{x, y) = S{x — y) with G{x, y) = 0 for y G 914). 

Let X be the n-algebra generated by the random variables v{x)(j)i for f G {1,..., m} 
(with V solution of (3.1)). We will identify the interpolation basis elements by condition¬ 
ing the solution of (3.1) on X. Observe that the covariance matrix of the measurement 
vector v{x)4>i)i<i<rn is the m x m symmetric matrix 0 defined by 


4'i{x)G{x,y)4>j{y) dxdy (3.3) 

Note that for I G I'^Ql = where w is the solution of (1.1) with right hand 

side g = Therefore (since the test functions 4)i are linearly independent) 0 

is positive definite and we will write 0“^ its inverse. Write 6ij the Kronecker’s delta 
((5i,i = 1 and 6ij = 0 for z 7 ^ j). 

Theorem 3.1. Let v be the solution of (3.1). It holds true that 

m 

E[z;(x)|T'] = '^'ijjiix) 
i=l 

where the functions tfi G are defined by 


v{y)(t>i{y) dy 


(3.4) 


'ifi{x) :=E v{x) / v{y)(l)j{y)dy = 6ij, j e {1,. ..,m} 
Jn 

and admit the following representation formula 


Mx) = / G{x,y)4>j{y)dy . 

i=i 


(3.5) 

(3.6) 


Furthermore, the distribution of v conditioned on T is that of a Gaussian field with mean 
(3.4) and covariance function T{x,y) = G{x,y) -|- YlTj=i 
- E™ 1 In dz - Ya=i In G(x, z}c/)i(z) dz . 
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Proof. The proof is similar to that of [77, Thm. 3.5]. The identification of the covari¬ 
ance function follows from the expansion of T{x,y) = E (u(x) — E[u(x)|T']) (u(y) — 

E[t;(y)|T']) . Note that (3.6) proves that ifi G Hq{Q). □ 

Since, according to (3.5) and the discussion preceding (3.1), each is an elementary 
gamble (bet) on value of the solution of (1.1) given the information = 6ij for 

j = 1,... ,m we will refer to the basis functions as gamblets. According to 

(3.4), once gamblets have been identified, they form a basis for betting on the value of 
the solution of (1.1) given the measurements (/f^ 

3.2 Optimal recovery properties 

Although gamblets admit the representation formula (3.6), we will not use it for their 
practical (numerical) computation. Instead we will work with variational properties 
inherited from the conditioning of the Gaussian field v. To guide our intuition, note that 
since C is the precision function (inverse of the covariance function) of v, the conditional 
expectation of v can be identified by minimizing given measurements constraints. 

This observation motivates us to consider, for i G {1,... ,m}, the following quadratic 
optimization problem 


(3.7) 


{ Minimize ||'0||a 

Subject to V' £ -f^d(^) = dij for j = 1 ,..., m 

where HV’IU is the energy norm of '0 defined in (1.3). 

The following theorem shows that (3.7) can be used to identify ifi and that gamblets 
are characterized by optimal (variational) recovery properties. 

Theorem 3.2. It holds true that (1) The optimization problem (3.7) admits a unique 
minimizer il^i defined by (3.5) and (3.6) (2) For w G M™, I'he unique 

minimizer of \\fi\\a subject to f^'tf{x)(l)j{x) = wj for j G {l,...,m} and (3) (using the 
scalar product defined in (3.2)J = ©["/• 

Proof. Let w G and fiw = with fii defined as in (3.6). The definition 

of 0 implies that f^fiwix)4>j{x) = Wj for j G {l,...,m}. Furthermore we obtain by 
integration by parts that for all ip G Hq{Q), {fiw, = YlTj=i In Therefore, 

if V’ G is such that il){x)(fj{x) = Wj for j G {1,..., m} then — 'f’w) = 0 

and 

Ufa = UwWl + U - TpwWl ( 3 - 8 ) 

which finishes the proof of optimality of fii and fiw D 


3.3 Optimal accuracy of the recovery 

De fin e 

m „ 

u*{x) :='^fii{x) / u{y)(j)i{y) dy 


(3.9) 


2 = 1 
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where u is the solution of (1.1) and V’i are the gamblets defined by (3.5) and (3.6). 
Note u* corresponds to Player IPs bet on the value of u given the measurements 
dy)i<i<rn- In particular, if v is the solution of (3.1) then 

u*(x) = E[u(x)| f v{y)(j)i{y) dy = [ u{y)(pi{y) dy] (3.10) 

Jq J u 

For (j) G write C~^(j) the solution of (1.1) with g = 4>. The following Theorem 

shows that u* is the best approximation (in energy norm) of u in span{£“^(/>j : i G 

{l,...,m}}. 

Theorem 3.3. Let u he the solution of (1.1), u* defined in (3.9) and (3.10). It holds 
true that 

||n-u*||a= inf \\u-fi\\a (3.11) 

'ip£spa,n{C~^ (pi:iG{l,... ,m}} 

Proof. By Theorem 3.1 span{£~^(/)j : i G m}} = span{xjji,... ,'iprn} and (3.11) 

follows from the fact that f^iu — u*)(f)j = 0 for all j implies that u — u* is orthogonal to 
spanl-f/;!,..., firn} with respect to the scalar product □ 


3.4 Transfer property and selection of the measurement functions 

We will now select the measurement (test) functions fii by extending the result of Propo¬ 
sition 2.10 to the continuous case. For V, a finite dimensional linear subspace of II~^(fl), 
define 

(divaV)“^y := span{(div oV)“^(/> : (3.12) 

where (divaV)“^()) is the solution of (1.1) with g = —(f). Similarly define A~^V := 
span{A“^(/> : </> G V} where A~^(fi is the solution of (2.10) with g = —(f). 

Proposition 3.4. Ifu andu' are the solutions of (1.1) and (2.10) (with the same r.h.s. 
g) and V is a finite dimensional linear subspace of then 


— , inf u —u < inf 

\/Amax(a) J^eA-W “ •ue(divaV)-W 


\u—v 


la < 


\/~Ki 


- inf 

fa) aeA-iy 


u —v\ 


(3.13) 


Proof. Write G the Green’s function of (1.1) and G* the Green’s function of (2.10). 
Observe that for / G P and u = (div aV)“^/, ||u — u||q = J^ 2 {g{x) — f{x))G{x,y){g{y)— 
f{y)) dxdy. The monotonicity of Green’s function as a quadratic form (see for instance 
[12, Lemma4.13]), implies f^ 2 {g{x)-f{x))G{x,y){g{y)-f{y))dxdy < f^ 2 {g{x)- 

f{x))G*{x, y){g{y) — /(y)) dx dy (with a similar inequality on the l.h.s.) which concludes 
the proof. □ 


This extension, which is also directly related to the transfer property of the flux-norm 
(introduced in [14] and generalized in [103], see also [115]), allows us to select accurate 
finite dimensional bases for the approximation of the solution space of (1.1). 
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Construction 3.5. Let (rj)i<i<m a partition of ri such that each ti is Lipschitz, 
convex and of diameter at most H. Let {4>i)i<i<m be elements of L‘^{Ll) such that for 
each i, the support of (pi is contained in the closure of Ti and (pi ^ 0. 

Proposition 3.6. Let (<^i)i<i<m be the elements of Construction 3.5 and let u be the 
solution of (1.1). If V = span{(/)j : 1 <i < m} then 


inf < CH\\g\\L2^n) 

t;G(div aV)~^V 


(3.14) 


with C = (vr.y/Amii((a)) (1 + maxi<j<m ) (writing |rj| the volume of Ti). 


Proof. Using Proposition 3.4 it is sufficient to complete proof when a is the constant iden¬ 
tity matrix. Let u' be the solution of (2.10) and v G A~^V. Note that Av = Cifpi, 
therefore \\u' - = - /^(m' - v){g - Ya=i Taking a = gf (pi we ob¬ 
tain that j^_{g — Cj(pj) = 0 and, writing |ri| the volume of Ti, \\u' — = 

“ Ya=i ^ M ~ '^))(5 “ which by Poincare’s inequality (see 

[91] for the optimal constant I/tt used here) lead to 11'^^ “^ ^ Y(aLi ( fr |V(tt' — 

n)P) ^ ( It id ~ ^ ■ Therefore, by using Cauchy-Schwartz inequality and sim¬ 

plifying, \\u' - ?^||h 1 (q) < fWd - YT=i Ci </’*!! 1,2(0) • Now, since each (pi has support in n 


we have || Y.T=i = T.T=iiU 

concludes the proof. 


< \\9 


|2 

Il 2 (o) 


<(,2 

maxi<i<m which 




□ 


The value of the constant C in Proposition 3.6 motivates us to modify Construction 
3.5 as follows. 


Construction 3.7. Let {4>i)i<i<m be the elements constructed in 3.5 under the addi¬ 
tional assumptions that (a) each (pi is equal to one on Ti and zero elsewhere (b) there 
exists 6 G (0,1) such that for each i G {1, ..., m}, Ti contains a ball of diameter 6H. 

Let {(pi)i<i<m be as in Construction 3.7. Note that the additional assumption (a) 
implies that the constant C in Proposition 3.6 is equal to 2/(tt. Assumption 
(b) will be used for localization purposes in subsections 3.5 and 3.6 (and is not required 
for Theorem 3.8). The following theorem is a direct consequence of Proposition 3.6 and 
Theorem 3.3. 


Theorem 3.8. If u is the solution of (1.1) and {ipi)^i 
(3.5), (3.7) and (3.6) then 

inf lilt —u||a< - -j^=^= 

Despan{V'i,...,'i/'m} 7 ryAmin(a) 

and the minimum in the l.h.s of (3.15) is achieved for v = 


are the gamblets identified in 

^I|5'IIl2(o) (3.15) 

u* defined in (3.9) and (3.10). 
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Remark 3.9. The assumption of convexity of the subdomains Ti is only used to derive 
sharper constants via Poincare’s inequality for convex domains (without it, approxima¬ 
tion error bounds remain valid after multiplication by tt). Similarly, the transfer property 
can be used to derive constructions that are distinct from 3.5 and 3. 7. 

Remark 3.10. Gamblets defined via the constrained energy minimization problems (3.7) 
are analogous to the energy minimizing bases of [65, Ilf, 122, 121] and in particular 
[108]. However they form a different set of basis functions when global constraints are 
added: the (total) energy minimizing bases of [65, Ilf, 122, 121, 108] are defined by 
minimizing the total energy W'fiWt subject to the constraint ^ related to 

the local preservation of constants. Numerical experiments [122] suggest that total en¬ 
ergy minimizing basis functions could lead to a 0{^/H) convergence rate (with rough 
coefficients). Note that (3.7) is also analogous to the constrained minimization problems 
associated with Polyharmonic Splines [5f, 30, 31, 32, 87], which can be recovered with a 
Gaussian prior (on f^) with covariance function 5{x — y) (corresponding to exciting (3.1) 
with white noise). We suspect that the basis functions obtained in the orthogonal decom¬ 
position of [6f] can also be recovered via the variational formulation (3.7) by identifying 
the null space of the Glement quasi-interpolation operator with that of appropriately cho¬ 
sen measurement functions 4>i. 

3.5 Exponential decay of gamblets 

Theorems 3.2 and 3.3 show that the gamblets ifi have optimal recovery properties anal¬ 
ogous to the discrete case of Theorem 2.2 and Corollary 2.6. However one may wonder 
why one should compute these gamblets rather than the elements (divaV)~^())j since 
they span the same linear space (by the representation formula (3.6)). The answer lies 
in the fact that each gamblet ifi decays exponentially as a function of the distance from 
the support of 4 >i and its computation can therefore be localized to a subdomain of 
diameter 0(77In ^) without impacting the order of accuracy (3.15). Consider the con¬ 
struction 3.7. Let ifi be defined as in Theorem 3.2 and let Xi be an element of r^. Write 
B{x,r) the ball of center x and radius r. 

Theorem 3.11. Exponential decay of the basis elements fji. It holds true that 
[ aN/'ifi < e^~T^ [ {VififaV'ifi (3.16) 

with 1 = 1-1 (e/7r)y^Amax(a)/Amin(a)(l + 22 (2/(5)^^'^/^) (where e is Euler’s number). 

Proof. Let k,l G N* and i G {1,..., m}. Let Sq be the union of all the domains tj that 
are contained in the closure of B{xi, klH) n H, let S'! be the union of all the domains Tj 
that are contained in the closure of {B{xi, {k -|- l)lH)f n H and let S'* = 5 q H S'f n H (be 
the union of the remaining elements tj not contained in Sq or Si). Let rj be the function 
on H defined by t]{x) = dist(x, So)/(dist(x, So) -|-dist(x, Si)). Observe that (1) 0 < ry < 1 
(2) rj is equal to zero on So (3) rj is equal to one on Si (4) || Vr/ll^oo^Q) < j^. Observe that 
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-/f^77V^idiv(aVV'i) = f^V(r]'ipi)a-i/ji = J^r](y^i)'^aV'ipi + Therefore 

fgJy'ipi)'^aVil;i < h + h with 

Ii = \\Vv\\L^^n){Yl [ ^i)H [ W*)^aVV^OVAmax(a) (3.17) 

and I 2 = — J^r]il!idiv{a'V'tpi). By (3.6), — div(aV'0i) is piecewise constant and equal to 
on Tj. By the constraints of (3.7) J)., ipi = 0 for i 7^ j. Therefore (writing rjj the 
volume average of rj over Tj) we have 

^2 < - ^ f {r]-r]j)'^|Jidiv{aV^pi) <j ^ ||V^i||i2(^^)|| div(aVV’i)llL2(^.). (3.18) 

TjCSlUS* TjCS* 

We will now need the following lemma 
Lemma 3.12. If v £ span{'0i,... ,V’m} then 

II div(aVu)||i2(f,) < H-^\\v\\a{XmM2^+^/6^+^)'2 (3.19) 

Proof. Let c G and v = YllLi CjV’i- Observing that — div(aVu) = YllLi ™ 7?' 

and using the decomposition || div(aVu)||^2(Q) = YllLi II div(aVu)||^2(.,-,), we obtain that 

m m 

||div(aVu)||22(^) = J^(J^Ci0j)Vil (3-20) 

i=l i=l 

Furthermore, v can be decomposed over Tj as u = ui +U2, where vi solves — div(aVui) = 
1 in Tj with ui = 0 on dTj, and V 2 solves — div(aVu2) = 0 in Tj with 

U2 = u on dTj. Using the notation |^|^ = observe that |Vu|^ = |Vui|q + 

f l^^2|a' Furthermore, f , |Vui|q = f vi. Writing Gj the Green’s func- 

'J '3 ' j 

tion of the operator — div(aV-) with Dirichlet boundary condition on dTj, note that 
vi = U©)"/) /r2 ^j(^> y) dx dy. Using the monotonicity of the Green’s function 

as a quadratic form (as in the proof of Proposition 3.4), we have f 2 Gj{x,y) dx dy > 

’’’j 

y— f ^2 G*{x,y) dx dy where G*j is the Green’s function of the operator —A with 
Dirichlet boundary condition on dxj. Recall that 2 J(_ G*j{x,y) dy is the mean exit time 
(from Tj) of a Brownian motion started from x and the mean exit time of a Brownian 
motion started from x to exit a ball of center x and radius r is (see for instance 
[12]). Since Tj contains a ball of diameter 6H, it follows that 2 f 2 G*{x,y) dx dy > 

{SH/4:)‘^~^'^Vd (where Vd is the volume of the d-dimensional unit ball). Therefore (after 
using |rj| < Vd{H/2)'^ and simplification), 

/ m 

|Vui|2 > (^c,0ri)V,|i72d2+V(2'+"A^ax(a)), (3.21) 

J i=i 

which finishes the proof after taking the sum over j. □ 
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Now observe that since ipi = 0 for i ^ j, we obtain, using Poincare’s inequal¬ 
ity (with the optimal constant of [91]), that ||V’i|lL2(rj) < \\'^'4’i\\L2{Tj)H/TT. Therefore, 
combining (3.17), (3.18) and the result of Lemma 3.12, we obtain after simplification 


f < ^VAmax(a)/Amin(a)(l + 2i(2/,5)i+''/2^ (3.22) 

J S\ J s* 


Taking I > ^ Amax(a)/Amin(o)(1 + 25 (2/(5)^++^) and enlarging the integration domain 
on the right hand side we obtain aV'ipi < e~^ fs*uSi We conclude 

the proof via straightforward iteration on k. □ 


3.6 Localization of the basis elements 


Theorem 3.11 allows us to localize the construction of basis elements 'ipi as follows. For 
r > 0 let Sr be the union of the subdomains tj intersecting B{xi,r) (recall that Xi is an 
element of r*) and let be the minimizer of the following quadratic problem 


Minimize (VaVip 

Subject to Ip £ HQ{Sr) and fg cpj'i’ = ^i,j for J such that Tj C Sr- 


(3.23) 


We will naturally identify with its extension to Hq (fl) by setting = 0 outside 
of Sr- From now on, to simplify the expression of constants, we will assume without loss 
of generality that the domain is rescaled so that diam(f2) < 1. 


Theorem 3.13. It holds true that 


W'tpi - < Ce 2m , (3.24) 

where I is defined in Theorem 3.11, C = (Amax(a)/\/\nin(^)F7~2“22^'’*+®/(\/T^(j'^+^) 
and Vd is the volume of the d-dimensional unit ball. 

Proof. We will need the following lemma. 

Lemma 3.14. It holds true that 

UiWa < (775)-5-VAmax(a)25''+2(Prf)-5 (3.25) 

where Vd is the volume of the d-dimensional unit ball, and, 

< e-3^77-2-'^A^ax(a)2™/(Pd5''+2) (3.26) 

where I is the constant of Theorem 3.11 and Vij is the distance between Ti and Tj. 

Proof. Since Ti contains a ball B{xi,6H/2) of center Xi £ Ti and diameter <577/2, there 
exists a piece-wise differentiable function r/, equal to 1 on B{xi,5H/4), equal to 0 on 
{B{xi, 5H/2)Y and such that 0 < ?? < 1 with || Vr/||j;^oo(Q) < Since ip = 
satisfies the constrains of the minimization problem (3.7) we have llV'ijla + llV’lla) which 


19 







proves (3.25). Theorem 3.2 implies that Observing that — div(aVV’i) 

is piecewise constant and equal to 0“^?^ on Tj and applying (3.21) (with v = and using 
|Vui|2 < |Vu|2), we obtain that 

(3.27) 

which leads to (3.26) by the exponential decay obtained in Theorem 3.11 and (3.25). □ 


0-/I < (A„,ax(a)25+'^/(52+^|r,|))^77-i( / 


Let us now prove Theorem 3.13. Let 5*0 be the union of the subdomains Tj not 
contained in Sr and let Si be the union of the subdomains Tj that are at distance at 
least H from Sq (for 5o = 0 the proof is trivial, so we may assume that S'o 7^ 0, similarly 
it is no restriction to assume that Si 7^ 0). Let r] be the function on 11 defined by 
ri{x) = dist(x, 5'o)/(dist(x, 5o) + dist(x,5i)). Observe that is a piecewise differentiable 
function on 11 such that (1) rj is equal to one on Si and zero on Sq (2) ||Vry||2,tx>(f2) < ^ 
and (3) 0 < r/ < 1. Since satisfies the constraints of (3.7), we have from (3.8), 


/ loc.r II 2 

L 


^loc,r ||2 


kWl 


(3.28) 


Let V'fc*' be the minimizer of [V^Y'aVilj subject to V’ £ ^oi^r) and fg (^jip = 5k,j 
for Tj C Sr- Write Wj = Let ■= Noting that = 

+ E., Wjil^Y 1 where S* is the union of Tj C Sr not contained in Si, and using 
property (3) of Theorem 3.2 (with 0^’^/^ = fg follows that 




/ loc,r ||2 


+ 


E 

r,CS* 


WjY 


+ 2 


E 

FiCS* 


0 




W-i 


(3.29) 


Noting that rjipi £ 77g(5r), Theorem 3.2 implies that ||V’E||a < IlhV'ilU) which, combined 
with (3.29) and (3.28) leads to \\ipi - “ llV’illa “ 2Er,cS* 

(using WrjipiWl - \\ipi\\l < fs* {vi’iYaV 

llV'i - V’|°''’nia < [ Vir]'i/;iYaV{r]il>i) + 2\ ^ ■ (3.30) 

Now observe that A Jg^V{r]'ipiYaV{r]'il;i) < Js* 

Applying Poincare’s inequality we obtain fg*\i^Y ^ Er^cS*/r (since 

= 0 for Tj C S*), and Jg, \A\^ < Ln(i?(x„r-2i7))-(Wi)'^aVV’i. Com¬ 

bining these equations with the exponential decay of Theorem 3.11 we deduce 

j V(?7'0i)^aV(7/V’j) < 2(1 + Amax(a)/(7r^Amin(a)))e^"^^||V’i||E (3-31) 

Similarly, using Cauchy-Schwartz and Poincare inequalities we have for Tj C S*, 

\wj\ < iTjl^WipiWL^irj) < I'Ll ^ (fr^ (VV’i)'^a(Vt^i))^/y/Amin(o) and 
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^ I Er,c5*(0ij- )^killM/5*(W*)'^a(VV’i)/Amin(a))2. Using (3.27) 


weobtainthat I - (Vax(a)2^+'^/(5^+'^) ^77 ^(/_g*(VV’-’’’)'^aVV'-’'’)' 

which by the exponential decay of Theorem (3.11) (and WtpiWa < llV’i’^IU) leads to 




(3.32) 


Using (3.25) to bound llV'l’^lla and combining (3.32) with (3.31) and (3.30) concludes 
the proof. □ 


The following theorem shows that gamblets preserve the 0{H) rate of convergence (in 
energy norm) after localization to sub-domains of size 0(771n(l/77)). They can therefore 
be used as localized basis functions in numerical homogenization [5, 85, 64, 87]. Section 
4 will show that they can also be computed hierarchically at near linear complexity. 

Theorem 3.15. Let u he the solution of (1.1) and (V’i°'^’^)i<i<m the localized gamblets 
identified in (3.23), then for r > H{Ci In ^ C 2 ) we have 

inf ll'W —i’||a< — / • (3.33) 

^ r /loc,r / loc,ri . Tfi\ 

The constants are Ci = {d + 4:)l and C 2 = 2/In 22^+^ ^ where I is the constant of 

Theorem 3.13. Furthermore, the inequality (3.33) is achieved for v = 

Proof. Let vi := nnd V 2 = YlT=i with Ci = ufii. Theorem 3.8 implies 

that ||n-ui||a < 2/{Tr^/Xfi^ffi(aj)H\\g\\L2(^Qy Observe that ||n - U2||a < ||w - uilja-b 

ll^^i — V 2 \\a and jjui — V 2 \\a < maxj \\fii — |ci|- Using Poincare’s inequality 

11^11^2(0) < diam(fl)||Vn||j;^2(Q) (with diam(ll) < 1) we obtain k*l — /n I'^l — 

1 

ll5llL2(n)2“‘^/^U^^ /Amin(o)- We conclude using Theorem 3.13 to bound maxj jjV’i— 

□ 


4 Multiresolution operator decomposition 

Building on the analysis of Section 3, we will now gamble on the approximation of the 
solution of (1.1) based on measurements performed at different levels of resolution. The 
resulting hierarchical (and nested) games will then be used to derive a multiresolution 
decomposition of (1.1) (orthogonal across subscales) and a near-linear complexity mul¬ 
tiresolution algorithm with a priori error bounds. 

4.1 Hierarchy of nested measurement functions 

In order to define the hierarchy of games we will first define a hierarchy of nested mea¬ 
surement functions. 
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Definition 4.1. We say that X is an index tree of depth q if it is a finite set of q- 
tuples of the form i = {ii,...,iq) with 1 < n < mo and I < ij < for 

j > 2, where mg and are strictly positive integers. For I < k < q and 

i = {ii,...,iq) G I, we write i^^^ := {ii,...,ik) and X^^'^ := : i G X}. For 

k < k' < q and j = {ji,... ,jk') G X^^ ) we write := (ji,... ,jk)- For i G X^^^ and 
k < k' < q we write i^^’^ ^ the set of elements j G X^^ ^ such that = i. 

Construction 4.2. LetX be an index tree of depth q. Let 6 G (0,1) and 0 < Hq < • • • < 
Hi <1. Let {Tf"\k G {1,... ,q},i G he a collection of subsets of Ll such that (1) 

for 1 < k < q, i G X^^^) is a partition ofLl such that each is a Lipschitz, convex 
subset ofLl of diameter at most Hk and contains a ball of diameter (2) the sequence 
of partitions is nested, i.e. for k ^ {1,... ,q — 1] and i G X^^\ := Ujgj(fc,fc+i) 

(k) 

As in Remark 3.9, the assumption of convexity of the subdomains ^ is not necessary 
to the results presented here and is only used to derive sharper/simpler constants. Let 

be the indicator function of the set (i.e. = 1 if x G and = 0 

(k)\ 

if X 0 r/ '). Note that the nesting of the domain decomposition implies that of the 
measurement functions, i.e. for k G {1,..., g — 1} and i G X^^\ 

^ (4.1) 

jgl(fc + l) 

where jg X^^'i x matrix defined by = 1 if j G and 

_^{k,k+i) _ Q j 0 ^(fc,fc+i)_ assume without loss of generality that = 

is constant in i (for the general case, rescale/renormalize each and the entries of 
7 r^k,k+i) corresponding multiplicative factors, we will keep track of the dependence 

of some of the constants on maxjj |T^^^^|/|rj^^|). 

4.2 Hierarchy of nested gamblets and multiresolution approximations 

Let us now consider the problem of recovering the solution of (1.1) based on the nested 
measurements (/^^for k G {1,... ,g}. As in Section 3 we are lead to inves¬ 
tigate the mixed strategy (for Player II) expressed by replacing the source term g with 
a centered Gaussian field with covariance function X = — div(aV). Under that mixed 
strategy. Player IPs bet on the value of the solution of (1.1), given the measurements 
Un'^iv)4>f\y) dy)i(ixi.>^), is (see Subsection 3.3) 

^xW(x) := y] V’f^(x) f u{y)4>f\y)dy, (4.2) 

iex('=) 

where (see Theorem 3.2), for k G {1,..., g} and i G X^^\ is the minimizer of 
f Minimize WifWa 

[^Subject to if £ ^ 13 (f^) and = dij for j G . 
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Define := Hq{Q) and, for /c G {1,..., g}, 

:= span{V^f ^ | i G (4.4) 

By Theorem 3.1 span{';/)|^^ | i G = span{£“^0j-^^ | i G and the nesting 

(4.1) of the measurement functions implies the nesting of the spaces The following 

theorem is (which is a direct application of theorems 3.3 and 3.8) shows that is the 
best (energy norm) approximation of the solution of (1.1) in 

Theorem 4.3. It holds true that (1) for k G q}, C and = 

span{£“^(/)|^^ I i G and (2) If u is the solution of (1.1) and defined in (4.2) 

then 

= inf ||M-u||a <'— 1 ^ ~ -^fc||g||L 2 (o) (4-5) 

vmw vr^Amin(a) 


4.3 Nested games and martingale/multiresolution decomposition 

As in Section 3 we consider the mixed strategy (for Player II) expressed by replacing 
the source term g with a centered Gaussian field with covariance function C. Under 
this mixed strategy, Player IPs bet (4.2) on the value of the solution of (1.1), given 
the measurements {f^u{y)(l)[^\y) dy)-^ 2 -(k), can also be obtained by conditioning the 
solution V of the SPDE (3.1) (see (3.10)), i.e. 

u(^)(x) = E u(a:) [ v{y)4>f\y) dy = [ u{y)(j)[’"\y) dy, i € (4.6) 

J J 

(k) 

Furthermore, each gamblet V’) ' represents Player IPs bet on the value of the solution of 
(1.1) given the measurements u{y)4>^^\y) dy = 6ij, i.e. 


= e[u [ v{y)(j)f\y) dy = 6ij, j G 


(4.7) 


Now consider the nesting of non-cooperative games where Player I chooses g in (1.1) 
and Player II is shown the measurements (f^ ^ ®tep by step, in a hierarchical 

manner, from coarse {k = 1) to fine {k = q) and must, at each step k of the game, gamble 
on the value of solution u. The following theorem and (4.6) show that the resulting 
sequence of approximations form the realization of a martingale with independent 
increments. 

Theorem 4 . 4 . LetFk be. the a-algebra generated by the random variables (f^ v{x)(l)'^'^)^^xW 
and 

u(^^(x) := E[u(x)|X'fc] = [ 'e{y)(t>f\y) dy (4.8) 

iex(fe) 

It holds true that (1) Fi,... ,Fq forms a filtration, i.e. Fk C Fk+i (2) For x G It, v^^\x) 
is a martingale with respect to the filtration {Fk)k>i, i-e. v^^\x) = E[u*^^'''^)(x)|Xfc] (3) 
and the increments — v^^'^)k>i are independent Gaussian fields. 
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Proof. The nesting (4.1) of the measurement functions implies C J-k+i and {J^k)k>i 
is therefore filtration. The fact that is a martingale follows from = E[u|j-fc]. 
Since and the increments — v^^^)k>i are Gaussian fields belonging to the 

same Gaussian space their independence is equivalent to zero covariance, which follows 
from the martingale property, i.e. for A: > 1 = E — 

=Epi)E[(u(^+^)-u(*=))|T'fc]j =0andforA;> j > - 

=e[(u(^'+^) -u(^))|T'fc]j =0. □ 


Remark 4.5. Theorem 4-4 enables the application of classical results concerning martin¬ 
gales to the numerical analysis of (and ). In particular (1) Martingale (concen¬ 
tration) inequalities can he used to control the fluctuations of (2) Optimal stopping 
times can he used to derive optimal strategies for stopping numerical simulations based 
on loss functions mixing computation costs with the cost of imperfect decisions (3) Tak- 
ing q = oo in the construction of the basis elements 'if] ^ (with a sequence Hk decreasing 
towards 0) and using the martingale convergence theorem imply that, for all (p G 
—>■ f^vip as k ^ oo (a.s. and in L^). 

The independence of the increments is related to the following orthogonal 

multiresolution decomposition of the operator (1.1). For defined as in (4.4) and for 
k G {2,..., g + 1} let be the orthogonal complement of within with 

respect to the scalar product (•, •)^. Write ®a the orthogonal direct sum with respect to 
the scalar product (•, •)^. Note that by Theorem 4.3, defined by (4.2) is the finite 
element solution of (1.1) in (in particular we will write = u). 


Theorem 4.6. It holds true that (1) For k G {2,..., g + 1}, 

Qj{fc) ^ Qj(i) 0^ 22 j(2) 0^ ... 0^ , (4.9) 


(2) for A: G {1,..., g}, belongs to and 

u = + • • • + (n*''^^ — n*''^“^^) + (u — u*-'^^) (4-10) 


is the orthogonal decomposition of u in (12) = ©a • • • ©a ©a , 

and (3) is the finite element solution of (1.1) in 

Proof. Observe that since the are nested (Theorem 4.3) belongs to 

Furthermore (by Property (1) of Theorem 4.3 and integration by parts), for 
i G belongs to span{ \ i G Finally, 

(4.2), the constraints of (4.3) and the nesting property (4.1) imply that for i G 

= Y)jei(k+i) 0 which implies that 

^^(^+1) _ y{k) belongs to 2B^^'''^\ □ 
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4.4 Interpolation and restriction matrices/operators 

Since the spaces 53^^^ are nested there exists a x matrix such that 

for 1 < A: < g — 1 and i G X^^) 

jgl(fc + l) 


We will refer to as the restriction matrix and to its transpose ;= 

{Rik,k+l))T 

as the interpolation/prolongation matrix. The following theorem shows that 
(see Figure 1) is Player IPs best bet on the value of given the infor¬ 
mation that = Si^s, s G 

Theorem 4.7. It holds true that for i G and j G , 

j^{k,k+i) ^ ^{k)^(k+i) = ]E[ v{y)(j)f^^\y) dy\ v(y)(/)j^\y) dy = 6^, I G X(^)]. 

Proof. The first equality is obtained by integrating (4.11) against and using the 

constraints satisfied by in (4-3). For the second equality, observe that since Xfe 

is a filtration we can replace v in the representation formula (4.7) by (as defined by 
the r.h.s. of (4.8)) and obtain 

which corresponds to (4.11). □ 

4.5 Nested computation of the interpolation and stiffness matrices 

Let V be the solution of (3.1). Observe that (f^ is a Gaussian vector with 

(symmetric, positive definite) covariance matrix 0^^^ defined by for i,j G X^^\ 

&ij -■= 4'f'\x)G{x, y)(ff\y) dxdy. (4.12) 

As in (3.3), is invertible and we write Q^k), i j^g inverse. Observe that, as in 

(k) 

Theorem 3.2, -0) ' admits the following representation formula 

^f\x)= X] / G{x,y)f)f \y)dy (4.13) 

jexW 

Observe that, as in Theorem 3.2, where A^^'> is the (symmetric, 

positive definite) stiffness matrix of the elements i.e., for i,j G 

Write ■nl-k+'^'k) transpose of the matrix T:^k,k+^) (defined below (4.1)) and 
the X(^) X X(^) identity matrix. The following theorem enables the hierarchical/nested 
computation of A^^'^ from 
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Theorem 4.8. For b G ig thg (unique) minimizer c G of 

{ Minimize 

Subject to Tr(^’'^+^)c = b 

FuTthcTUlOTC ^ ^ 

Q(k) ^ ^(k,k+l)Q(k+l)^(k + l,k) 

J^(k) _ j^{k,k+l) j^{k+l) j^{k+l,k) _ (4-16) 

Proof. Using the decompositions (4.11) and (4.1) in Igg^^Jg 

j^(k,k+i)^(k+i,k) _ j(fc)_ Using (4.13) and (4.1) to expand Theorem 4.7 leads to 
j^(k,k+i) _ ^(fc)^(fc,fc+i) 0 {fc+i)^ Using (4.1) to expand and in (4.12) leads to 
0 (^) = Using (4.11) to expand and in (4.14) leads to 

(4.16). Let b G Theorem 3.2 implies that Xliexf'') is the unique minimizer 

of ||u||^ subject to u G Hq{Q,) and = bj for j G Since C 

and since the minimizer is in the minimization over v G iLg(fl) can be reduced to 

V G of the form v = X]iei('=+i) which after using (4.1) to expand the 

constraint = bj, corresponds to (4.15). □ 

4.6 Multiresolution gamblets 

The interpolation and restriction operators are sufficient to derive a multigrid method 
for solving (1.1). To design a multiresolution algorithm we need to continue the analysis 
and identify basis functions for the subspaces For k = 2,... ,q let be the 

finite set of /c-tuples of the form i = (ii ,... ,ik) with 1 < ii < mo, 1 < ij < 
for 2 < j < k — 1 and 1 < ifc < — 1. Where the integers m. are the same 

as those defining the index tree X. For a matrix M write Im(M) and Ker(M) its image 
and kernel. 

Lemma 4.9. For k = 2,.. . ,q let be a x X^^^ matrix such that Im(lU^*^^’^) = 
Ker(7r(*^“^’*^)). It holds true that the elements iXi^'^)i^j(k) £ defined as 

xf'’ ■■= E (4.17) 

jei(fe) 

form a basis 

Proof. Since = span{(/iE^^ I i £ w G 23^^^ belongs to if and only if 

ffi = 0 for all j G which, taking w = and using (4.1), translates into 

Writing \ the number of elements of (which is equal to 
the dimension of observe that Therefore Im(lU(^)’^) = 

Ker( 7 r(*^“^’*^)) also implies that the \ elements linearly independent and, 

therefore, form a basis of 2 B^^^. □ 
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Remark 4.10. Observe that since 0 = it also 

holds true that and = Ker(R(*^-i’*^)) . 


From now on we choose, for each k G {2,... ,g}, a x matrix as in 

Lemma 4.9. This choice is not unique and to enable fast multiplication by (or its 
transpose) we require that for (j,i) G xZ^^\ = 0 if Therefore, 

the construction of requires, for each s G to specify a number — 1 of 


ms-dimensional vectors y ..., ^ that are linearly independent and 

orthogonal to the ms-dimensional vector (1,1,...,1,1). We propose two simple con¬ 
structions. 


(k) 




Construction 4.11. For k G {2, ...,g}, choose such (1) = 0 for {j,i) G 

j'W X with ^ and (2) for s G t G {!,... ,rns — 1} and t' G 

{ 1 , • • • ,'ms}, 

For k € {2,..., q} and i = {h,..., ik-iAk) G define z+ := (ii,..., ik-iAk + 1) 
and observe that under construction 4.11, 

xf ^ = V'f ^ - V’S^ (4.18) 

whose game-theoretic interpretation is provided in Figure 1. 


(a) 


T,' 




W T 


(k) 
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Figure 1: If {ts^\ s G is a nested rectangular partition of fl then (a) is Player 

IPs best bet on the value of the solution u of (1.1) given f (k) u = 6ij for j G X^^^ (b) 

b 

is Player IPs best bet on u given f (k) u = 6ij — Si+ j for j G X^^^ (c) ig 

b ’ ’ 

Player IPs best bet on /^(fe+i) u given f^(k) u = 6ij for j G X^'^) . 


For the second construction we need the following lemma whose proof is trivial. 

Lemma 4.12. Let be the sequence of n x n matrices defined (1) for n = 2 by 

= (1,-1) and = (1,1) and (2) iteratively for n > 2 by for 
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1 < *, j < JT', ^n +1 j' ~ 1 1 < j < n + 1 , = 0 /or 1 < i < re — 1 arec? 

= —re. Then for re > 2, 7/ie rows of are orthogonal, = 1 for 1 < j < n 
and we write i/ie corresponding orthonormal matrix obtained by renormalizing the 
rows of . 

Construction 4.13. For k G {2,...,q}, choose such (1) = 0 for {j,i) G 

jF) X with and (2) for s G and t G {l,...,rres — 1} and 

t' G {1,.. . ,ms} (where fs defined in Lemma f.l2). 

Observe that under Construction 4.13 (1) the complexity of constructing is 

\j^(k-i)\ ^ ^2 ^ 2 ) nr(Opjr(OX = j{k) -^^ere is the x identity matrix. 


4.7 Multiresolution operator inversion 

( 1 ) ik) 

We will now use the basis functions if I ' and y • to perform the multiresolution inversion 
of (1.1). Let be the x (stiffness)matrix observe 

that 

(4-19) 

Observe that is positive, symmetric, definite and write B^^'l’ ^ its inverse. Let 
^(fc,A:+i) 2:(fc) ^ j(A:+i) j^^atrix defined by 




^(/c,fc+l) 


/(^(fc,fc+i)^{fc+i,fc)).^. 


(4.20) 


Using the notations of Definition 4.1 note that = nii. Let 

be the x matrix defined as 


jj(k,k-l) ._ _Q{k)-ly^{k)j^(k)-{k,k-l) 


and write := its transpose. 

Theorem 4.14. It holds true that for /c G {1,..., g — 1} and i G 



,{k) _(fc,fc+l) , (fc+1) p,(fc,fc+l) (fc+1) 

V’i = 2^ K,i + 2^ DP >x) ’■ 

;gl(fe+i) 

(4.22) 

In particular, 

j^{k,k+l) _ -(fc,fc+l) _|_ ^(A;,fc+1)^(A:+1) 

(4.23) 

Proof. For s G X*^^) write := X)«el('=+i) and := 

X(^)}. Let X G y G and 

span{'i/)i^^ s G 


if='^ Xsiff^ + ^ ■ 

(4.24) 


5gi(fc) jeX(''+i) 
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If = 0 then integrating V’ against for i G (and observing that = 

6i^s) implies x = 0 and y = 0. Therefore the elements form a basis for 

Observing that dim(QJ*^^'''^^) = dim(Tl(^)) + dim(2n^^'''^^) we deduce that 
^(fc+i) _ ^(k) _j_ Therefore, since C can be decomposed as in 

(4.24). The constraints Jq = ^i,s lead to Xg = 5i,s. The orthogonality between 

V’ and leads to the equations (V’) = 0 for j £ i.e. 

)a + = 0 ) '^'focli foans- 

+ that is (4.22). Plugging (4.17) in (4.22) and 

comparing with (4.11) leads to (4.23). □ 

Let g be the r.h.s of (1.1). For k G {1,..., g} let be the |Xl^l|-dimensional vector 
defined by gf^^ = J^i’i^^'g for i G X^^l. Observe that g^^^ can be computed iteratively 
using 

g{k) ^ j^{k,k+i)g{k+i) _ (4 25) 

For k G {2 ,... ,q}, let be 1771*^11-dimensional vector defined as the solution of 

^(fc)y;(0 = (4.26) 

Furthermore let t/1^1 be the |Xl^l|-dimensional vector defined as the solution of 

= 5 ( 1 ) (4.27) 

According to following theorem, which is a direct consequence of Theorem 4.6, the 
solution of ( 1 . 1 ) can be computed at any scale by solving the decoupled linear systems 
(4.26) and (4.27). 

Theorem 4.15. For k G {2,... ,q}, let be the finite element solution of (1.1) in 
It holds true that = Yli&jW in particular, 

+ E E (4.28) 

igX(l) fc'=2 jgj-(fe') 

4.8 Uniformly bounded condition numbers across subscales/subbands 

ik) 

Taking g = 00 in Theorem 4.6, the construction of the basis elements ipi ' leads to the 
multiresolution orthogonal decomposition. 

Hi (54) = ©a 2rr(*). (4.29) 

i=2 

In that sense the basis elements and x'^'^ could be seen as a generalization of wavelets 
to the orthogonal decomposition of Hl{Xl) (rather than L^(n)) adapted to the solution 
space of the PDF (1.1). We will now show that this orthogonal decomposition induces 


(fc+i) 


V _(/C,/C + l)/ y (/C+l) 

lates into 
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a subscale decomposition of the operator — div(aV) into layered subbands of increasing 
frequencies. Moreover the condition number of the operator — div(oV) restricted to 
each subspace will be shown to be uniformly bounded if Hk-i/Hk is uniformly 

bounded (e.g. if is a geometric sequence). Write Hq := 1 and let 5 be defined as in 
Construction 4.2. 


Theorem 4.16. If k £ {1,... ,q} and v £ then 

:l+d/2 

:Hk < 


\/Amax(a)25/2+'^/2^^'' ' || div(aVu) ||i 2 (n) ’ 

If {k = 1 and v £ or {k £ {2,, q} and v £ then 


< 


div(aVu)||i2(f^) yX 

min 


--Hk- 


1 • 


(4.30) 


(4.31) 


Proof. (4.30) is a direct consequence of Lemma 3.12. For k 
consequence of Poincare’s inequality. Let k £ {2,...,g}. 
Theorem 4.3 imply 


|p||a ^ f 

SUp^g 2 IJ('') II div(aVj;) 11^2(0) - ™D'e27('=-i) || div(aV?;)|L 2 (n) 


< 


1 (4.31) is a simple 
©„ 2 rr(^) and 


's/ .^min (^) 


Hk-1. 


□ 


Write |c| the Euclidean norm of c and for A; G {1,..., g} let 


71 . := inf 




and 7 fc := sup 




ceM- 


i(fe) 


(4.32) 


- I fA:)i I (fc) I 

Write |t| the volume of a set r and note that 7 ^ < max^gj(fc) |ri | and 7 ^ > min-gj(fe) |rl ^ |, 
therefore ^kl'Jk < 

For a given matrix M, write Cond(M) := ^^/M)/\/ A min (M'^M) its condi¬ 
tion number. 


Theorem 4.17. It holds true that 


Cond(Xl(^)) < 


1 An,ax(a)25+^ 
Hi An,in(a),52+2d’ 


and for k £ {2,... ,q}, 


Cond(i?W) < (%i)2(^Lp4^)2|^Cond(W('=)w('=)’^). 
Hk Amin (a) 


(4.33) 


(4.34) 


Furthermore, Cond{W^^^W^^^’'^) = 1 under Construction 4-13 andCon(l{W^^^W^^^''^) < 
2(^Hk-i/{5Hk))^^ under Construction 4-11. 
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Proof. Let k G {!,...,(?} and c G Write v = ■ Observing that 

ll^lla = c^A^'^'>c and || div(aVn)||2 2 ^^^ = || E*eiW ^> jk\A^>^^c\'^, (4.30) 

implies that /{Xraa.x{o,)2^~^^) < c^A^^'>c/\A^^'>c\‘^, which, after taking the mini- 

mnm in c leads to (for k > 1) 

Amax(^(")) < An,ax(a)25+V(^|<5'+S), (4-35) 


and for A: > 2 (using (4.19)) 

Amax(i?^'=^) < A„,ax(kL('=)w('=)’'^)A„,ax(a)25+V(^|<5'+S) • (4-36) 


Similarly for A; = 1 (4.31) leads to Amin(^*'^^) > Amin(o)/ 7 i- Now let us consider 
k G {2,... ,q} and c G Write w = (4-17) and (4.19) 

imply that ||tc||„ = and (using (4.32)) 

||div(aVu;)b 2 (^) = ||Eiej(.),,- 6 xw(AlWwW’^c),<))f Ob¬ 
serving that w G (4.31) implies that )pV(fc)>^c - A~R^fc-i- Taking 

c = for y G RT'") deduce that < % x^)Hl_v Writing 

.^g have obtained that 

Amin(a)/(iL|_i7feAmax(iV('^)Tjv(fc))) < . (4.37) 

For yt G {2,..., q} let := T^{k,k-i) ji{k-i,k) _ Ug^^g ^(fc-LO = ^(fc-i)7r(fc-bO0(O and 
^(fc-i,fc)0(fc).^(fc,fe-i) _ 0(fc-i) ^Theorem 4.8) we obtain that (p(0)2 = pl^)^ i.e. p(0 ig 
a projection. Write ||fW ||Ker( 7 r(fe-i.'=)) := sup,j,gKer(P'=-i.'')) 

Lemma 4.18. It holds true that for k G {2,... , q}, 


Amax(A^^")’^iV('^^) < 


1+ l|T(^)|lKer(^(fe-i,fc)) 

Amin(LF(fc)W(fc).^) 


(4.38) 


Proof. Since Im(kF^^)Tj and Im(7r^^A i)) are orthogonal and dim(M^^*’^) = dim(Im(lF^^^Tj)_|_ 
dim( Im( 7 r(^A-i) ^ fQj- X G there exists a unique y G and z G such 

that X = _|_|^|2 _ |p|/(A:),ry |2 _|_ |,^(fc,fc-i)^| 2 ^ Observe that = 

W^k)wik),Ty (since VF(*^)7r(^A-i) ^ pl^-hO^; = ^(fc-i,fc)^(fc),ry ^ ^ 

j^(fc-i,fc).^(fc,fc-i) _ j(fc-i) fj-Qm Theorem 4.8). Therefore, |xp = _|_ |p(fc)(x — 

iT(OTy ^|2 .^ii^]^ y — (yY(k)yY(k),T'^-iy^(k)^^ ^ Taking x = and 

observing that P^^^x = 0 (since Jl(k i,k)j^{k)]Y{k),T _ g fj-gm the (•,•)^-orthogonality 
between and 211^^^) leads to |— ^■^/{k),Ty^ 2 _^^p(k)'^{k),Ty ^2 .^it]^ y _ 

(WWwWT)-ip(O^. Therefore |^WwWT^;| 2 < (i+||p(0 

which concludes the proof after taking v = and maximizing the l.h.s. over 

\v'\ = 1. □ 
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(4.39) 


Lemma 4.19. Writing ||M ||2 := supx\Mx\/x the spectral norm, we have 


lipwii 


2 


< 


xGKer(7r(^“^’^)) 


x'^Q(^')x 

x'^x 


Proof. Let x G Ker( 7 r(^“^’^)). Using we obtain that 

\p(k)x\ = || 2 |( 0 (*^)) 2 x|. Observing that for 

M = {qW )^2 we have MM^ = 0(*'“^) and for N = 2 

we have Xmaxi^'^N) = XmaxiNN'^) we deduce 

||.^(fc,fc-i)^(fc-i)^(A:-i,fc)^ 0 (A:)j 21|2 _ ||.jj-(A:,fc-i)^(A:-i).^(fc-i,fc)||^ conclude by taking the 
supremum over x G Ker( 7 r(^“^’^)). □ 


Lemma 4.20. It holds true that 


sup 


^Qik) 


X ^^9 

“ ^ ^k-l 


fk 


yfc7r^Amin(®) 


(4.40) 


Proof. Let y G and a G M. Let x = Write (f = YliexW Xi4>\^^ and 

^|J = (—div(aV-))“^()). Observe that ||V’||a = x'^&^^^x > Us¬ 
ing /q = 0 for z / / and selecting a = ||)) 2 (q) (assuming, without loss 

of generality, that |li 2 (Q) = is constant in i, for the general case, rescale 


each by a multiplicative constant) we obtain that for j G J( 


d^-i) - 


n 

2 


EiexW bOvL(^)’^y)j = 0. Therefore, since W'lfWl = 

we have for if' G span{(/.f"^^ | i gX^^-^)} ||V^||2 = < UW L^n)U - W l^q) ■ 

Choosing if' = Eiex('=-i)Ili 2 (f 2 ) obtain (via Poincare and 
Cauchy-Schwartz inequalities as in the proof of Proposition 3.6) that \\if — if^Wi^^n) < 
Hk-i\\if\\a/ (vT-yAmin(a)) and deduce ||V’||a < Hk-i\\4>\\L2(^Q)/{7r^/XfnwJaf)- Observing 
that ||<('||| 2 (f 2 ) — and 7 ^ < a~^ < 7 ^ we summarize and obtain that 

< ^ 2 _J^(fc),Ty| 2 ^ 2 y'('^^^ 2 _)^^.^('^^)^ which 

concludes the proof of the lemma (since Ker( 7 r^^“^’^^) = Im(lU^^^’'^)). □ 

Observing that < A ma y(7r(^’^~^)7r^^~^’^))A ma. ^(^*-^~^^) and 

using (4.35), we derive from lemmas 4.19 and 4.20 that 

7i25+‘^A„ 


||p(0|| 


Ker(7r('“-i''=)) — '^max(7r' ’ 'TT' 


r(A:A-l)„(fc-l, 


fc)^ 


_ c(«) 

7fc7fc-l(52+'^7r2Amin(a) ’ 


(4.41) 


Observing that jg block-diagonal and using the notations of Definition 4.1 

we have Amax(7r(''’^“^)7r(^“bfc)) = max^g 2 -(fc_i) sup^.gjj’"^ I = max^g 2 -(i,-i) mj. 

Noting that a set can contain at most (maXjg 2 :(fc_i) |rj^ ^^|)/(min^g 2 :(fc) sub¬ 

sets rf^'^ we have 

max mj < [Hk-i/{5Hk)Y (4.42) 
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and conclude that Amax(7r^^’^ !.*:))< (^Hk-i/{6Hk))'^. Therefore (4.37) and Lemma 

4.18 imply, after simplification, that 




Amin(o) 






(4.43) 


Recalling that jkHk < using jkllk-i < Bp{Hk-i5Y, and summarizing we con¬ 
clude the proof of (4.33) and (4.34). Recall that under construction 4.13 we have 
VL(^)VL(^)T = j(fc) which implies Cond(lT(^^lT(^^’'^) = 1. Under construction 4.11, 
iy(^)ll^(^)T ig block diagonal with for j G diagonal blocks corresponding to 

{vfij — 1) X (nij — 1) matrices such that (1) for n = 1 and x G M, x^M^^'>x = 2x^ 

(2) for n = 2 and x G M^, = x\ + (x 2 — xi)^ -|- x| and (3) for for n > 3, 

and X G M"', x^M*^”)x = x\ + YPi=iPi ~ Note that for all n > 1, 

Amax(-A7(”')) < 3. Furthermore, for n > 3 (n < 2 is trivial), introducing the variables 
2/2 = X 2 - Xi,... ,yn = Xn - Xn -1 we obtain that x'^M^^'>x = x^ + y 2 + ■■■ +Vn + x"^ 
and |xp = xf -f (xi -b 2 / 2 )^ + ■ ■ ■ + (xi + y 2 + ■ ■ ■ + ynY < (x'^M("')x)n(n -b l)/2. 
Therefore, Amin(-Tf^”^) > 2l{n{n -b 1)). We conclude that under construction 4.11 
Cond(iy('')iy('=)’'^) < max^gj{fc_i) 3(mj — l)mj/2 and bound rrij as in (4.42). □ 


4.9 Well conditioned relaxation across subscales 

If Hk is a geometric sequence or if Hk-i/Hk is uniformly bounded, then, by Theorem 
(4.17), the linear systems ((4.26) and (4.27)) entering in the calculation of the gamblets 
(and therefore and the subband/subscale solutions and 

have uniformly bounded condition numbers (in particular, these condition numbers are 
bounded independently from mesh size/resolution and the regularity of a(x)). Therefore 
these systems can be solved efficiently using iterative methods. One such methods is the 
Conjugate Gradient (CG) method [57]. Recall [98] that the application of the CG method 
to a linear system Ax = h (where A\s a, n x n symmetric positive definite matrix) with 
initial guess yields a sequence of approximations x^^^ satisfying (writing |e|^ := 

e^^e) |x - x(')|a < ~ where Cond(^) := Amax(^)/Amin(^)- 

Recall [98] also that the maximum number of iterations required to reduce the error by 
a factor e (|x — x^^'^\a < (■\x — x*^*^)!^) is bounded by 2y^Gond(4.) In | and has complexity 
(number of required arithmetic operations) 0(-^Gond(A)A^^) (writing Na the number 
of non-zero entries of A). 


4.10 Hierarchical localization and error propagation across scales 

Although the multi-resolution decomposition presented in this section leads to well con¬ 
ditioned linear systems, the resulting matrices B^^'> and A^^'> are dense and to achieve 
near-linear complexity in the resolution of (1.1) these matrices must be truncated by 
localizing the computation of the basis functions (and therefore The ap¬ 

proximation error induced by these localization/truncation steps is controlled by the 
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exponential decay of gamblets and the uniform bound on the condition numbers of the 
matrices To make this control explicit and derive a bound the size of the localiza¬ 

tion sub-domains we need to quantify the propagation of truncation/localization errors 
across scales and this is the purpose of this subsection. 

For k G g}, p > 1 and i G we define ( 1 ) as the subset of indices 

j G whose corresponding subdomains are at distance at most H^p from 
and (2) S'* := Let pi,...,pq > 1. For k G — 1} and i G Z^^\ 

write as the subset of indices j G such that G i^. For i G let 

^(g+i):loc Pq], k £ {1,... ,q} and i G let be the minimizer of 

Minimize Ill'll a subject io ij} £ ^^d f = 5ij for j £ (4.44) 

Jn 

where for A; G g — 1} and i £ Z^^\ ig defined (via induction) by 

Tjf+i)’l°c := span{^Af| j G 

From now on we will assume that for some H £ (0, 1) (or simply that 

is uniformly bounded from below and above by H^). To simplify the presentation, we 
will also, from now on, write C any constant that depends only d, fi, Ainin(fl), A mav la). 6 
(e.g., 2C'Amax(a) will still be written C). The following theorem allows us to control the 
localization error propagation across scales. For k £ {1,..., g}, let pg ^ 

X(^) matrix defined by := and let £{k) be the (localization) 

error £{k) := ( EiexW llV'f ^ “ V’f " • 

Theorem 4.21. For k £ {l,...,g-l}, we have 

£{k) < CH-i£{k + 1 ) + Ce"^'=/^i/5-(fc+b(rf+i), 45 ^ 

Proof. We will need the following lemma summarizing and simplifying some results 
obtained in Theorem 4.16 when ■ 

Lemma 4.22. Let and be as in Construction f.ll or Construction 4.13. 

It holds true that fork G {g,...,2} (1) ||iy (^)||2 < ^3 (2) l/Amin(VF(^)VF(^)’^) < CR-^^^ 

(3) ||7r(^-^’'')||2 < CRi (4) ||7r(*=-iA)||2 < CR-^/^ (5) II 2 < CR-^/^ (6) 

Cond{B^^y) < CR-^-^<^ (1) Amax(5^*^^) < 

(8) 1/Amin(5('')) < Ci/^( 2 +d)- 2 - 2 d_ Furthermore, (9) Cond(^(i)) < CR-^ (10) 1/Amin(^(^^) < 
CR'^ andfork£{l,...,q} (11) < CR-^^^+^\ 

Proof. From the proof of Theorem 4.17 we have (1) and l/Amin(kF^^^kF^^^’^) < 
max-g 2 :(fe_i)(mj — l)mi/2, which implies (2). For (3), noting that = 0 if ^ i 

and 7 f|^- = l/rui otherwise, we have II 2 = max^g 2 ;{fc_i) < CR^. (4) 

follows from A max (7r*^^’^~^^7r^^~^A)j = max^g 2 :{*;-i) m, < CR Let us now prove (5). 

Using (4.23), (4.21) and defining as in Lemma 4.18, we have 

^ ^{k-l,k) ^ ^(k-l,k) ^{k)^^{k) ^ Ig^^g ^ ^ 
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Iblllb)- Using Lemma 4.18 and (4.41) we obtain that < 

(1 + and therefore < CR-'^. 

Summarizing we have obtained (5). ( 6 ), (7), ( 8 ) and (11) follow from Theorem 4.16 
and in particular (4.36), (4.43) and (4.35). See (4.33) and the proof of Theorem 4.17 for 
(9) and (10). □ 


We will also need the following lemma. 

Lemma 4.23. Let k G {1,..., (7 — 1} and let R he the X matrix defined by 

and Rij = 0 for j G j ^ ^(k+i) jj^p^,k+i ^ holds true that 

||i?||2 < CR^I‘^e-P'‘l^. 


Proof. Observe that \\R\\l < maxX)jgi(fc+i)/iPfe,'=+i 1-^1 with < 

Let i G Using Theorem 4.7 and Cauchy-Schwartz inequality we have 

|^(fc,fc+i)| ^ ||V;f^||^ 2 (^jfe+i))ll<bf^^^lli 2 (.,)fe+i))- Therefore ^^.g 2 ;(fc+i)/iPfc,fc+i I< 
CH^+iEj&x(^+^)/iPk.k+^ \\^l’'^\\l 2 ^fik+i)y Observe that Ej^x(^+^)/iPk,k+^ ^ 

Jls&xW/iPk = 0 for s / i we obtain from Poincare’s in- 

equality that ^ Therefore ^ 

Using Theorem 3.11 we obtain that 
J2s&xW/iPk ^ C'e-^'^/’blb’f^lla- Using (3.25) we have 

therefore X)j 6 x('=+i)/iPfc’'=+i I- T'77'^e“‘^ □ 


Let us now prove Theorem 4.21. We obtain by induction (using the constraints 
in (4.44)) that for k G {l,...,g} and i G X^^b satisfies the constraints of 

(4.3). Moreover (3.8) implies that if fi satisfies the constraints of (4.3) then 


-( fc )||2 


+ 


Ha- 


Therefore, for A: G {2, ...,g — 1}, is also the mini- 

mizer of Wfi — over functions if of the form fi = Y2j^iPk<f‘+'^ Pj ^(fc+l),l0C gg^l^jgfyij^g 

the constraints of (4.44). Thus, writing fi* := Ylj&iPk^k+i we have 

(since ip* satisfies the constraints of (4.44)) < IIV’* “ Write 

V’l := EjeKfe+i) and ^^2 := YljexC^+i)/iPk,k+i Observ¬ 
ing that b* = V’i“V ’2 we deduce that lla ^ 2 ||' 0 i—V’i^^||a+ 2 ||b’ 2 ||a after 

summing over i, {£{k))‘^ < 2 (/i -h h) with R = Eiei('=) II Ejei(''+i) - 

V’f’^i)’^°c)||a and h = EiexW II Ejei('=+i)/i^fc’'“+i Writing S the 

2 ;(fc+i) xX(*’+^) symmetric positive matrix with entries Sij = ^ “ 

note that R = Trace[i?(^T+i) 572 (fc+i.*:)j, Writing the matrix square root 
of 5, observe that for a matrix U, using the cyclic property of the trace, Trace[17517^] = 
Trace[ 52 C/'^f 7 S' 2 ] < Amax(U’^U) Trace[5], which (observing that Trace[5] = (X(A: -|- 1))^ 
and Amax(U^U) = \\U\\l) implies R < ||i?*^^’^'''^bl 2 (T(^ + 1))^- Therefore (using Lemma 
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4.22) we have ^/li < CH 2 £{k + 1). Let us now bound l 2 - Let R be defined as 
in Lemma 4.23. Noting that j^g^yg (-gg g^ove) 

I 2 = Trace[i?j4(^+^)’*°‘^i?^] < Amax(-R^.R) Trace[^(*^'''^)’*°'^]. Summarizing and using 
Lemma 4.23 we deduce that £{k) < CH 2 £{k + 1) + CH 2 e“^*/^-\/Trace[^(^+^)’*°‘^]. Ob¬ 
serving that Trace [^(A:+i),iocj < _j_ _j_ y^Trace[A(^+^)] and using Trace[^^^^^^] < 

maXjg 2 :{fc+i) Ila (Lemma 3.14) lU < C'-f^fc+i ^5 we conclude the 

proof of the theorem. □ 

Let be the finite element solution of (1.1) in := span{'0j^^’^°‘^ | j G 

For k G {2,..., g}, let be defined as in Construction 4.11 or Construction 

4.13. For i G let := EjeK^) For fc G {2,... ,g} let - 

.g(fc-i),ioc i^g ^]^g fjgj^g element solution of (1.1) in := span{Xj^^’'°'^ | j G 

For A: G {2,... ,g}, write ^(i),loc - u^^-^k^ocy Let fifOloc 

the X matrix defined by := . Observe that = 

W(k)Amocw(k),T_ wLte for k G {2,...,q}, £{k,x) := (E.ejW llxf - xf"• 
The following theorem allows us to control the effect of the localization error on the 
approximation of the solution of ( 1 . 1 ). 

Theorem 4.24. It holds true that for k ^ {2,... ,q\ (1) £{k,x) H CH~'^/‘^£{k). Fur¬ 
thermore for /c G {2,..., g} and £{k, x) < C~^we have 
(2) Cond(S('=)’'°‘=) < C'iT- 2 - 2 '^, and (3) - (^Wloc _ ^(fc-i),ioc)||^ < 

C£{k,x)\\g\\H-^{n)H'^^^^‘^^‘^^-^‘^~^. Similarly for £{1) < we have 

U) Cond(AF),ioc) < CH-'^, and (5) ||■u(F _ < C£{l)\\g\\H-nn)H-‘^^<^/‘^■ 

Proof. We will need the following lemma. 

Lemma 4.25. Let xi, • • ■, Xm be linearly independent elements of Hq{£L). Let Xi, •. •, Xm 
he another set of linearly independent elements of Hq{Q). Write £ := (X^I^illXi “ 

x(lla)^- Fet B (resp. B') he the m x m matrix defined by Bij = {xiiXj)^ (resp. 
^ij — (Xi)Xj)(jJ- Fet Um (resp. u'^) be the solution of ( 1 . 1 ) in spanjxi | i = 
l,...,m} (resp. span{x( | i = l,...,m}). It holds true that for £ < y/\fahJj^/2 
(1) Cond(S') < 8 Cond(S) (2) \\B - B'h < 3y/Xrae^{B)£ (3) \\B-^ - {B')-% < 
12x'Xn...{B){Xn.iniB))-^£ and (4) \\um - u'JU < CX||5||^-i(o)^f^. 

Proof For (1) observe that x^KffffW) = sup|^|=i || X)™ 1 ^JiXilla < \/Amax(-B) -h £ 
and \/Xmi„(B') = vafx\=i\\YT=i^iX'i\\a > \/Amin(-B) - £. For (2) observe that for 
G M”* with \x\ = |y| = 1 we have y^{B - B')x = {YT=iyiiXi - Xi), Ei^i - 

( Ya=i Vix'i, E™ 1 XiiXi - Xi))a < (\/Amax(-B') -h Xras^{B))£. (3) follows from (2) and 
||S"^-(fi')"^ll 2 < ||S-.B'||2/(Amin(-B)Amin(.B'))- For (4) observe that Um = YnLi'^iXi 
(resp. u'^ = Ei^i ^iXi) where w = B~^h with bi = J^gXi (resp. w' = {B')~^b' with 
K = In9Xd- Therefore \\um-u'j^\\a < \w\£+ \w- w'\^/X((^J(B). w-w' = B-^{b-b')- 
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B — B')w' leads to |u) — w'\ < C'(||fl'||H-i(n)‘^ + II-® “ -B'|| 2 |'u^^|)/Amin(-B). Using (2), 

Amin(5)|u;|2 < WYJlLiWiXiWl < \W\\l < and Xmm{B')\w'\‘^ < 

we conclude the proof of (4) after simplification. □ 

Let us now prove Theorem 4.24. Using 
and noting that = 0 for ^ have (T(fe,x))^ < YhjejW 

Therefore, {£{k,x)f < 

(T(/c))^ max-g 2 ;{fc-i) ruj max^gjr(*,) Observing that (see (4.42)) 

max-g 2 -(fe_i) mi < 1/{H5Y and Zlieif'') - Amax(lU(*^)pU(^)’^) < 3 (see Lemma 

4.22) we conclude that (1) holds true with C = (3/5“^) 2 . (2) and (3) are a direct 

application of lemmas 4.25 and 4.22. For (3), observe that (resp. — 

u(^-i)’i°c) is the finite element solution of (1.1) in (resp. 211^^^’*°'^ := span{xj^^’^°'^ | 
j G The proof of (4) and (5) is similar to that of (2) and (3). □ 

Theorem 4.26. Let A; G {1,..., g}. lUe have, £{k) < C Yl'j=k . 

Proof. By Theorem 4.21, for fc G {1,..., (? — !}, we have £{k) < ak + bk£{k + l) with ak = 
(jg-pk/C and bk = CH~ 2 . Therefore we obtain by induction that £{k) < 

a-k + bkttk+i + bkbk+iak+2 H-h 6fc • • • bq- 2 aq-i + bk - ■■ bq-i£{q). Using Theorem 3.13 

we have £{q) < Ce~P‘i^^ and obtain the result after simplification. □ 

Theorem 4.27. Let eG (0,1). It holds true that if pk > (^((l + la ^ + In ^) 

for k £ {1,... ,q} then (1) for A: G {1,..., g} we have — u(^Aioc||^ < e||g'||j:^-i(Q) 
and ||u - < C{H^ + e)\\g\\L2(n) Cond(T(i)’i°") < CB-’^, and for k G 

{2,...,q} we have (3) Cond(i?*^^)’^°'^) < and (4) — 

Proof. Theorems 4.3 and 4.24 imply that the results of Theorem 4.27 hold true if for 
k G {!,...,(?} £{k) < (7-i/7-fc(i+'i/2)+7(i/2+3g^^2_ Ugiag Theorem 4.26 we deduce that 
the results of Theorem 4.27 hold true if for k G {l,...,g} and k < j < q we have 
Ce~^ 2 /^C^~^H~ 2 B 2 -p 2 < }{-kp+d/2)+7d/2+3^i^f^2j2y -^Ve conclude after simplihca- 
tion. □ 


5 The algorithm, its implementation and complexity 

5.1 The initialisation of the algorithm 

To describe the practical implementation of the algorithm we consider the (finite-element) 
discretized version of (1.1). Let 7/i be a regular fine mesh discretization of £l of resolution 
h with 0 < /i <C 1. Let M be the set of interior nodes Zi and N = |AA| be the number of 
interior nodes {N = 0{h~'^)) of Th- Write ((^OisA/" ^ set of regular nodal basis elements 
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(of ifg (O)) constructed from Th such that for each i G M, support((^j) C B{zi, CqK) and 
for y G 

jh%\‘^ < II < ih^\y\^ (5.1) 

i&M 

for some constants j,^,Co ^ 0{1). In addition to (5.1) the regularity of the finite 
elements is used to ensure the availability of the inverse Poincare inequality 

II^^IIl2(q) < Cih ^||u||i2(Q) (5.2) 

for V G span{(/?j | i G Af} and some constant Ci ^ 0(1), used to generalize the proof of 
Theorem 4.16 to the discrete case. 

Given g = Yli&M 9ifi want to find u G span{</?j | i G Af} such that for all j G Af, 

Wj^u) = / (fjg for all j € Af (5.3) 

Jn 

In practical applications a is naturally assumed to be piecewise constant over the fine 
mesh (e.g. of constant value in each triangle or square of Th) and one purpose of the 
algorithm is the fast resolution of the linear system (5.3) up to accuracy e G (0,1). 

u 

0.15 

0 1 

0.05 

0 

■0.05 


Figure 2: The (fine) mesh Th, a (in logj^g scale) and u. 




Example 5.1. IFe will illustrate the presentation of the algorithm with a numerical 
example in which Th is a square grid of mesh size h = (1 + 2 '?)“^ with q = 6 and 
64 X 64 interior nodes (Figure 2). a is piecewise constant on each square of Th and 

given by a{x) = OLi (l +0.5 cos (2^7r(^ + ^))^ (^1 + 0.5 sin (2^7r(^ - 3^))^ 

for X G ^ [ 29 ^’ contrast of a (i.e., when a is scalar, the ratio 

between its maximum and minimum value) is 1866. The finite-element discretization 
(5.3) is obtained using continuous nodal basis elements ipi spanned by {1, xi, X 2 , X 1 X 2 } in 
each square ofTh- Writing Zi the positions of the interior nodes ofTh, we choose, for our 
numerical example, g{x) = (cos(3ziq + 2 : 4 ^ 2 ) + sin(32;j^2) + sin(72;iq — f>Zi^2))Ti{‘^)- 
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Figure 3: and 


The first step of the proposed algorithm is the construction of the index tree X 
of Definition 4.2 describing the domain decomposition of Definition 4.2. To ensure 
a uniform bound on the condition numbers of the stiffness matrices (4.19) one must 
select the resolutions Hk to form a geometric sequence (or simply such that Hk-i/Hk 
is uniformly bounded), i.e. Hk = for some H G (0,1) (for our numerical example 
H = 1/2, q = 6 and we identify Z^^'> as the indices of the interior nodes of a square 
grid of resolution (1 + 2^)“^ as illustrated in Figure 3). In this construction H‘^ = h 
corresponds to the resolution of the fine mesh and each subset {i G X^'^^) contains 
one and only one element of J\f (interior node of the fine mesh). Using this one to one 
correspondence we use the elements of X = X'^ to (re)label the nodal elements 
as {(pi)i^x- The measurement functions are then identified (1) by selecting 

= ifi for i G X^'^) and (2) via the nested aggregation (4.1) of the nodal elements (as 
commonly done in AMG), i.e. for 

fe G {1,..., g — 1} and i G Z^^\ 

Remark 5.1. We refer to Figure 4 for an illustration of these measurement functions 

(k) 

for our numerical example. Note that the support of each c/i is only approximatively 
(and not exactly) and that the are only approximate set functions (and not 
exact ones). This does not affect the design, accuracy and localization of the algo¬ 
rithm presented here because the frame inequalities (4.32), and the Poincare inequalities 
WHi&iik) for x G Ker( 7 r(^“^’*=)), hold true. Indeed, 

(5.1) and Construction 4-2 imply that the frame inequalities (4.32) with 7 ^ < and 
Ik > and the Poincare inequalities are regularity/homogeneity conditions on the 
mesh and the aggregated elements. Although a fine mesh has been used to facilitate 
the presentation of the algorithm, the proposed method is meshless (it only requires the 
specification of the basis elements {ip/ji^x)- 
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Figure 4: The functions (j)^ with k G {1,..., g} and q = 6. 


5.2 Exact gamblet transform and multiresolution operator inversion 

The near-linear complexity of the proposed multi-resolution algorithm (Algorithm 2) is 
based on three properties (i) nesting (ii) uniformly bounded condition numbers (hi) local¬ 
ization/truncation based on exponential decay. Truncation/localization levels/subsets 
are, a priori, functions of the desired level of accuracy e G (0,1) in approximating the 
solution of (5.3) and to distinguish the implementation of localization/truncation (and 
its consequences) we will hrst describe this algorithm in its zero approximation error 
version (i.e. e = 0 and without using localization/truncation. Algorithm 1). Although 
this error-free version (Algorithm 1) performs the decomposition of the resolution of 
the linear system (5.3) (whose condition number is of the order of ^ 1) into the 

resolutions of a nesting of linear systems with uniformly bounded condition numbers, 
it is not of near linear complexity due to the presence of dense matrices. Algorithm 2 
achieves near-linear complexity by truncating/localizing the dense matrices appearing 
in Algorithm 1 (e-accuracy is ensured using the off-diagonal exponential decay of these 
dense matrices). Let us now describe Algorithm 1 in detail. Lines 1 and 2 correspond to 
the computation of the (sparse) mass and stiffness matrices of (5.3). Line 4 corresponds 
to the calculation of level q gamblets defined as the minimizer of ||V^||a subject to 
~ V’ £ span{<y9; | I G X}, note that since the number of constraints 

is equal to the number of degrees of freedom of ip, and since level q 

gamblets do not depend on a and are obtained by inverting the mass matrix in Line 3 
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Algorithm 1 Exact Gamblet transform/solve. 


1 

For i,j e Mij = (fiipj 

/ / Mass matrix 

2 

For i,j G Aij = J^(S/(pi)'^aVpj 

// Stiffness matrix 

3 

Compute M~^ 

/ / Mass matrix inversion 

4 

ForiGX(''), = 

/ / Level q gamblets 

5 

For i E = gi j j g^^^ 

= In 9 with g = Eiex('j) 9m 

6 

For i,j E X(5), aSJ = 

// ^(9) = M-^AM-^’’^ 

7 

for fe = (? to 2 do 


8 


U Eq. (4.19) 

9 

-uj(k) = ^(fc),-i]^(fc)^(fc) 

// Eq. (4.26) 

10 

For i E = EjelW 

// Eq. (4.17) 

11 

uW _ u(fc-i) = 

// Thm. 4.15 

12 

jjik,k-l) _ _J^ik)-l^/{k)^{k)^{k,k-l) 

// Eq. (4.21) 

13 

j^{k-l,k) — jj.{k-l,k) _|_ D(k-l,k)\Y{k) 

// Eq. (4.23) 

14 

j^{k-l) _ j^(k-Xk) ji^{k) j^{k,k-l) 

// Eq. (4.16) 

15 

For i E 

// Eq. (4.11) 

16 

g(k-l) _ j^{k-l,k)g{k) 

// Eq. (4.25) 

17 

end for 


18 

t/(i) = 

// Eq. (4.27) 

19 

= Eteim 

// Thm. 4.15 

20 

U = -|- -I- • • • -I- ^^) 

// Thm. 4.6 with u = 


(note that by (5.1), the mass matrix is of 0(1) condition number). Although not done 
here, one can also initialize the algorithm (and its fast version) with (which 

is equivalent to using level q measurement functions). Line 5 corre¬ 

sponds to initialization of the vector introduced above (4.25). Line 6 corresponds 
to the initialization of the stiffness matrix introduced in (4.14). The core of the 
algorithm is the nested computation performed (iteratively from k = q down to k = 2) 
in lines 8 to 16. Note that this nested computation takes and 

inputs and produces (1) and outputs for the next itera¬ 
tion and (2) the subband of the solution and subband gamblets {Xi^'^)i£j(k) 

(which, do not need to be explicitly computed/stored since Line 11 is equivalent to 

Note also that the gamblets (V'i*^^)jex('=) 

(Xi^'^)i^j(k) can be stored and displayed using the hierarchical structure (4.11). Through 
this section and the remaining part of the paper we assume that the matrices are 
obtained as in Construction 4.11 or 4.13. Note that the number of non-zero entries of 
^{k-i,k) ]Y{k) ig C)(|x(^)|) (proportional to H~^ in our numerical example). Lines 9 

corresponds to solving the well conditioned linear system = and the 

|X(^“i)| well conditioned linear systems Note that by 
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Theorem 4.17 the matrices have uniformly bounded condition numbers and these 
linear systems can be solved efficiently using iterative methods (such as the Conjugate 
Gradient method recalled in Subsection (4.9)). is computed in lines 19 and 20 (re¬ 
call that is also of uniformly bounded condition number) and the last step of the 
algorithm, is to obtain u via simple addition of the subband/subscale solution and 
(^(fc) _ 

~^^) 2 <k<q- Observe that the operating diagram of Algorithm 1 is not a V or 
W but an inverted pyramid (or a comb). More precisely, the basis functions are 
computed hierarchically from fine to coarse scales. Furthermore as soon as the elements 
have been computed, they can be applied (independently from the other scales) to 
the computation of (the projection of u onto corresponding to the 

bandwidth 



Figure 5: The basis elements ip^ with A: G {1,..., 6}. 


ik) 

Example 5.2. We refer to figures 5 and 1 for an illustration of the gamblets ip^ and 

corresponding to Example 5.1 with defined by Construction fill. We refer 

ik) 

to Figure 6 for an illustration of the exponential decay of the gamblets . ITe refer 
to Figure 8 for an illustration of the condition numbers of and (with W^^'l 
still defined by Construction f.ll). Observe that the bound on the condition numbers 
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Figure 6: Exponential Decay. 








Figure 7: The basis elements ipj and Xi with A; G {2,..., 6}. 


of depends on the contrast and the saturation of that bound occurs for smaller 
values of k under low contrast. We refer to Figure 9 for an illustration of the subband 
solutions ..., corresponding to Example 5.1. Observe that 

these (subband) solutions form a multiresolution decomposition ofu as a sum of functions 
characterising the behavior ofu at subscales [H'^, H],... Once the 

components — u^^'>,..., and have been computed one obtains, via 

simple summation, ..., the finite-element approximation ofu at resolutions 
H, , ..., illustrated in Figure 11. As described in Theorem f.3 the error of the 
approximation of u by u^^'i is proportional to for A: G {1,..., g — 1}. For k = q, as 


43 

















6 X 10^ ___ 

" An.ax(A('°)) 

3. An3ax(B<'^h 

A„3in(S(^)) 

2 ' 

,. High 
contrast 


Figure 8: Condition numbers of and 




Figure 9; — u^^\..., and 



illustrated in Figure 11, this approximation error drops down to zero because there is no 
gap between H'^ and the fine mesh (i.e., and tpi span the same linear space in the 
discrete case). Moreover, as illustrated in Figure 10, the representation of u in the basis 


formed by the functions 




( 1 ) 


m 


(i)i 


and 


Jt 

lixfil 


is sparse. Therefore, as illustrated in Figure 


11 one can compress u, in this basis, by setting the smallest coefficients to zero without 
loss in energy norm. 
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Figure 10: The coefficients of u in the expansion u 




+TI=2 


lixf'l 


5.3 Fast Gamblet transform/solve 

Algorithm 2 achieves near linear complexity (1) in approximating the solution of (5.3) 
to a given level of accuracy e and (2) in performing an approximate Gamblet transform 
(sufficient to achieve that level of accuracy). This fast algorithm is obtained by localiz- 
ing/truncating the linear systems corresponding to lines 3 and 12 in Algorithm 1. We 
define these localization/truncation steps as follows. For k G {l,...,g} and i G 
define as in Subsection 4.10 (i.e. as the subset of indices j G whose corresponding 
subdomains are at distance at most H]^p from r-^^^). 

Definition 5.2. For i G Z^'^\ let be the x F’' matrix defined by = Mij 

for l,j G . Let be the \iP‘>\-dimensional vector defined by = 6j^i for j G . 

Let y^'^’Pfi be the {iPi {-dimensional vector solution = gO-Pq). Wg define the 

solution of the localized linear system M’^'P'^M, = 5.^i (Line 3 of Algorithm 

2) as the -vector given by for j G 

Note that the associated gamblet (Line 4 of Algorithm 2) is also the solution 

of the problem of hnding fi G span{(/?j | j G iP‘>} such that firpj = dij for j G (i.e. 

localizing the computation of the gamblet to a subdomain of size HgPg). Line 5 
can be replaced by = g^ without loss of accuracy simplihes 
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Compression ratio = 105 
Energy norm relative error = 0.07 



Gamblet compression 


Figure 11: Relative approximation error in energy norm in log;^Q scale. 

Compression of u over the basis functions Xi^\ ■ • • i by setting 99% of the small¬ 
est coefficients to zero in the decomposition of Figure 10. 


the presentation of the analysis). Line 12 of Algorithm 2 is defined in a similar way as 
follows. 

Definition 5.3. Let k G {2,...,g} and B be the positive definite x matrix 
^(fc),ioc computed in Line 8 of Algorithm 2. For i G let p = Pk-i o,nd let i^ 

be the subset of indices j G such that G F (recall that if j = (R ,... ,ik) C 

then := (jd,..., Jfc-i) C be the i^ x matrix defined by 

= Bi j for l,j G i^. Let 6^’^) be the \i^\-dimensional vector defined by = 
— jor j G i^. Let yb.p) the \i'^\-dimensional vector solution 

of dFg define the solution of the localized linear system 

Inv(R(Gloc^(fe,fc-i),ioc ^ j{k) xx(fc-i) sparse matrix 

given by = 0 for j ^ and for j G i^. (Line 

13 of Algorithm 2) is then defined as the transpose of . 
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Algorithm 2 Fast Gamblet transform/solve. 


1 : For i,j e Mij = (fiipj 

2 : For i,j £ Aij = aVipj 


3: = 5. i 

‘ .L 1 






4: For i £ 

5: ForzEX('?),5l")’'“ = 4#’'°^5 
6: For i,j £ X(^), 
for /c = (? to 2 do 

Q(k),loc _ ^^(k)^{k),\oc-^{k),T 
y^(A:),loc _ ^^(fc),loc^-l^(fc)^(fe),loc 

For i G J("), = E,ei('=) 


// O(A') 
// O(A') 

// Def. 5.2, Thm. 5.6, C>(Apglnmax(i,g)) 


7 

8 
9: 

10 : 

11 : 

12 : 


,(fc),loc _ „(fc-l).loc ^ ^(fc),loc^(A:),loc 


// 

// 

// o(ivpf) 

// o(|x('=)|pf) 

// Thm. 5.6, C>(|X(*^)|pf Ini) 

// 0(|xW|pf) 

// 

// Def. 5.3, Thm. 5.6, 


13: 

14: 

15: 

16: 

17 

18 


pfc_i) 

0(|XW|ptiPfclni) 

^(fc-i,fc),ioc ^ ^(fc-i.fc) ^{fc-i,fc),ioc|^(fc) II 5 3^ C>(|X('=-i)|pf_^ 


^(/c—l),loc _ ^(fc—l,fc),loc^(fc),loc^(/c,/c—l),loc 

For i G X(^-i), = E,-6 x('=) 

^(fc—l),loc _ ^(/c—l,/i;),loc^(/i;),loc 

end for 

[j(l),loc _ ^(l),loc,-l^(l),loc 

19, = E.aw ' 

20: + (u(2)-l°= - u(l)iloc) + • • • + 


// o(|x(^-i)|pti) 
// o(|x(^-i)|pti) 

// 0(|X«|;9flni) 

// o(^pf) 
// 


Remark 5.4. Definition 5.3 (Line 8 of Algorithm 2) is equivalent to localizing the 
computation of each gamblet to a subdomain of size i.e., the gamblet 

fi^k computed in Line 15 of Algorithm 2 is the solution of (1) the problem of 

finding ip in the affine space X]jex('=) + span{xj'^^’^°'^ | G iP’^-^} 

such that Ip is (•, orthogonal to span{xj^^’^°'^ | £ iP*^-"^}, and (2) the problem of 

minimizing H'i/’IU in span{ip'^^^’^°^ \ l^^~^') £ iP>=-^} subject to constraints jQ4>f^ = ^i,j 
for j £ iP^-^. 

5.4 Complexity vs accuracy of Algorithm 2 and choice of the localiza¬ 
tion radii pk 

The sizes of the localization radii pk (and therefore the complexity of Algorithm 2) 
depend on whether Algorithm 2 is used as a pre-conditioner (as done with AMG) or 
as a direct solver. Although it is natural to expect the complexity of Algorithm 2 
to be significantly smaller if used as pre-conditioner (since pre-conditioning requires 
lower accuracy and therefore smaller localization radii) we will restrict our analysis 
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and presentation to using Algorithm 2 as a direct solver. Note that, when used as 
a direct solver, Algorithm 2 is parallel both in space (via localization) and in band- 
with/subscale (subscales can be computed independently from each other and 
and resolved in parallel). We will base our analysis on the 

results of Subsection 4.10 and in particular Theorem 4.27. Although obtained in a 
continuous setting, these results can be generalized to the discrete setting without dif¬ 
ficulty. Two small differences are worth mentioning. (1) In this discrete setting, an 
alternative approach for obtaining localization error bounds in the first step of the al¬ 
gorithm (the computation of the localized gamblets is to use the exponential 

decay property of the inverse of symmetric well-conditioned banded matrices [26] : since 
M is banded and of uniformly bounded condition number [26] (see also [11, Thm 4.10]) 
implies that decays like exp (— dist(r^^'^\ which guarantees that the bound 

£{q) < ( 7 //-'^/ 2 -'j( 2 +ci/ 2 )g-P 5 /C Theorem 4.26) remains valid in the discrete 

setting. ( 2 ) Since the basis functions (pi are not exact set functions, neither are the 

resulting aggregates This implies that, in the discrete setting, 

(k) 

necessarily equal to zero if rj ^ is adjacent to Sp^ (with j 0 , using the notation of 

Subsection 4.10). This, however does not prevent the generalization of the proof because 
the value of ^{k),ioc^{k) -g adjacent to S*^) can be controlled via the expo¬ 

nential decay of the basis functions (e.g. as done in the proof of Theorem 3.13). We will 
snmmarize this generalization in the following theorem (where the constant C depends 
on the constants Ci,Co,^ and 7 associated with the finite elements (ipi) in (5.1), in 
addition to d, Q, Aniin(a)j Amax(fl), <5). 

Theorem 5.5. Let u he the solution of the discrete system (5.3). Let — 

ri(^-i)’i°c, j^{k),ioc jj(k),\oc ^^g qJ Algorithm 2. Let and — 

^^(^-1) tiiQ outputs of Algorithm 1. For k G {2,...,g}, write ;= -|- 

^ ^ holds true that if pk > C((l-|- In ^-|- 

In A) for k e {I,... ,q} then (1) for /c G { 1 ,..., g — 1 } we have < 

e\\9\\H-An) and llu(^) - < C{H’^ + e)\\gh2(^n) (2) Cond(A(i)-i-) < CH-\ 

and for k G {2,... ,q} we have (3) Cond(i?*^^)’^°'^) < and (4) — 

^^(fc),ioc Finally, ( 6 ) 11 ^^ - < ell5ll^p-i(t^) . 

Therefore, according to Theorem 5.5 if the localization radii pk are chosen so that 
Pk = 0{^ lnmax(l/e, 1/Flk)) for k G {1, ... ,q} then the condition numbers of the matrices 

^(fc),loc ^(l),loc 

remain uniformly bounded and the algorithm achieves accuracy e in 
a direct solve. The following theorem shows that the linear systems appearing in lines 
3, 9 and 12 of Algorithm 2 do not need to be solved exactly and provide bounds on the 
accuracy requirements (to simplify notations, we will from now on drop the superscripts 
of the vectors y and b appearing in definitions 5.2 and 5.3). 

Theorem 5.6. The results of Theorem 5.5 remain true if (1) pk > C'((l+jj^pyj^) In |^-|- 
InA) for k G {l,...,(jr} (2) For each i G the localized linear system M^’^iy = 6.^i 
of Definition 5.2 and Line 3 of Algorithm 2 is solved up to accuracy \y — ^ 
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(using the notations of Subsection 4-9, i.e. |e|^ := e'^Ae, and writing 
y'^PP the approximation of y) (3) For k G {2,...,q} and each i G the localized 

linear system = b of Definition 5.3 and Line 12 of Algorithm 2 is solved up 

to accuracy \y — {k — 1)^. (4) For k G {2.,...,q} the 

linear system of Line 9 of Algorithm 2 is solved up to accuracy 

\y - y^pP|ij(fe).ioc < e||5||^^-i(n)/(2g). 


Proof. From the proof of Theorem (4.27) we need £{k) < C *:(i+<^/2)+7d/2+3g^^2 
for k G {1,... ,g}. By the inverse Poincare inequality (5.2) this inequality is satisfied 
ioT k = q for ^ for each i G which by 

the definition of M^’Pi and Line 4 of Algorithm 2 leads to (2). For k G {2, ...,g} 
the inequality £{k — 1) < — 1)^ is satisfied if for i G 

— 1)^. Using the notations of 

Definition 5.3 we have, + E,6«x 

with which leads to (3) by lines 15, 13, 10 and 8 of Algorithm 

2. For (4) we simply observe that for y G || Eie jw(y-rPP)*xf^’'°'^lla = \y- 

y"PP|BW,ioc. □ 


Let us now describe the complexity of Algorithm 2. This complexity depends on 
the desired accuracy e G (0,1). Lines 1 and 2 correspond to the computation of the 
(sparse) mass and stiffness matrices of (5.3). Note since A and M are sparse and 
banded (of bandwidth 2d = 4 in our numerical example) this computation is of 0{N) 
complexity. Line 3 corresponds to the resolution of the localized linear system introduced 
in Definition 5.2 using the x sub-matrix of M. According to Theorem 5.5, 

the accuracy of each solve must be \y — < C~^H'^^^'^~^^elq^. Since = 

0{pq) and since M®’P« is of condition number bounded by that of M, for each i the 
linear system of Line 3 can be solved efficiently (to accuracy using 

0{pq) = C>( lnmax(^, q)) iterations of the CG method (reminded in Subsection 4.9) with 
a cost of 0{pq) per iteration, which results in a total cost of C>(A^Pq lnmax(^, q)). Lines 

4 and 5 are naturally of complexity 0{Npg). Since = 0 if and are at a 

distance larger than 2H'^pq the complexity of Line 6 is O^Np^^^). Note that and 

^(fc),loc 

are banded and of bandwidth 0{Npf,). It follows that Line 8 is of complexity 
0(|X(^)|p^). According to Theorem 5.6 the linear system of Line 9 needs to be solved 
up to accuracy \y — y'^PPI^cfc) ,loc < e|blb-i(0)/2- Since is of uniformly bounded 

condition number this can be done using Ofln^) iterations of the CG method with a 
cost of C>(|X*^^)|p^) per iteration (using 0{\J^^'>\) = C>(|X(^)|)), which results in a total 
cost of C>(|X(^)|p^ In ^) for Line 9. Storing the fine mesh values of _ ^(fc-i),ioc 

in Line 11 costs 0{Npf) (since for each node x on the fine mesh only 0{pf.) localized 
basis functions contribute to the value of rt(^)4oc _ .y(A:-i),ioc^_ According to Theorem 
5.6, for each i G the linear system B^'^’P^y = b of Line 12 needs to be solved 

up to accuracy \y — < C~^H~^^'^^/‘^^‘^e/{k — 1)^. Since the matrix B^^’P'^ 

inherits the uniformly bounded condition number from this can be done using 
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0{ln^) iterations of the CG method with a cost of = 0{p'^_^pf.) per 

iteration. This results in a total cost of In for Line 12. We obtain, 

using the sparsity structures of j^{k-i,k),ioc complexity of Line 

13 is 0{\I^^-^^\pi_^H-'^) = and that of Line 14 is 0{\I^^-^^\pli-^pi). 

The complexity of lines 15 to 16 is summarized in the display of Algorithm 2 and a simple 
consequence of the sparsity structure of Line 18 is complexity 0{\I^^^\pf In 

(using CG as in Line 9). As in Line 11, storing the values of costs 0{Npf). Finally, 

obtaining in Line 20 costs 0{Nq) (observe that q = C)(lnA^)). 


Compute and store ^(fc),ioc^ ^(fc),ioc 

and s.t. ||u - ^ 49\\H-Hn) 

e < HI 

e > H‘1 

First solve 

Aln^'^i 

e 

Aln^-^ A 

Subsequence solves 

Aln“'+^ i 

e 

Aln"^ Alni 

e 

Subsequent solves to compute s.t. 

^ Ce\\g\\L2^^) 


Aln''+^ 1 

Subsequent solves to compute the coefficients 

of = J2iex(^) cS^Vf ^ 

s.t. \\u - < C'e(|| 5 ||i 2 (n) -h Ibllup) 


^—d 1 

e 

Subsequent solves to compute 

\\u - < Ce{\\gh2(n) + Iblkip) 


^~d — 

e 

a periodic/ergodic with mixing length < e, 

first solve of 

k “ < Ce{\\g\\]^2(^Q) + 5 Lip) 


(A(ln^‘^A)AP 

+e-4 

In^+i 1 
€ 


Table 1: Complexity of Algorithm 2. 


Total computational complexity, first solve. Summarizing we obtain that the 
complexity of Algorithm 2, i.e. the cost of computing the gamblets 
their stiffness matrices (^1^)’*°^^ and the approximation such that ||u — 

f0(A'( lnmax(^, A'd))^ ) (with Line 14 being the correspond¬ 
ing complexity bottleneck). The complexity of storing the gamblets (Xj^^’*°'^) 

and their stiffness matrices is 0(A( lnmax(^, . 

Computational complexity of subsequent solves with g € H ^(0). If (5.3) (i.e. 
(1.1)) needs to be solved for more than one g then the gamblets the 

stiffness matrices do not need to be recomputed. The cost of subsequent solves 

is therefore that of Line 9 i.e. C)(A( lnmax(f, Ad))^ In f) to achieve the approximation 
accuracy ||u - < e||ff||H-i(o)- 
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Computational complexity of subsequent solves with g G L?{Q) and e > 

If <7 G L‘^{Q) (i.e. if || 5 ||L 2 (r 2 ) is used to express the accuracy of the approximation) 
and e G [H^, H^~^] then, by Theorem 5.5, achieves the approximation accuracy 

||tt — < C'e||g||j;^ 2 (Q) (i.e. — itOTioc i^q j^g computed 

for j > k) and the corresponding complexity is 0(A^(ln (if g G then the 

energy of the solution can be in the fine scales and — liOi.ioc need to be 

computed for j > k). 


Computational complexity of subsequent solves with g Lipschitz continuous 
and e > H‘^. Note that the computational complexity bottleneck for computing the 
coefficients of in the basis ('0^^^’^°^) when g G C‘^{Q) and e G is in 

the computation of the vectors for j > k. li g is Lipschitz continuous then 

g{k)^oc approximated g{xf^^) where is any point in without loss of accuracy. 
Note that this approximation requires (only) evaluations of g and leads to a 

corresponding satisfying \\u - < Ce{\\g\\i 2 (^Q) + HffULip) (with H^llup = 

suPa;,j;go 19 ( 3 ;) “ 9{y)\/\^ ~ y\)- Therefore the computational complexity of subsequent 

solves to obtain the coefficients in the decomposition 

is 0(e“‘^(ln ^)'^'*“^) (i.e. independent from if o' is Lipschitz continuous). Of course, 
obtaining an iLQ(f2)-norm approximation of u with accuracy requires expressing the 
values of (and therefore ti(^)’*°'^) on the fine mesh, which leads to a total cost of 

0(A^(ln ^)'^). However if one is only interested expressing the values of on the fine 

mesh in a sub-domain of diameter e then the resulting complexity is 0((iVe'^+e“'^)(ln 


Computational complexity of subsequent L^-approximations with g Lipschitz 
continuous and e > Let (a;f)),gj(.) be points of (T|^^)jg 2 ;{'=) forming a regular 

coarse mesh of Q of resolution and write the corresponding (regular and coarse) 
piecewise linear nodal basis elements. If (as in classical homogenization or HMM) one is 
only interested in an L^-norm approximation of u with accuracy then the coefficients 

defined above are sufficient to obtain the approximation r' w 

in A 

that satisfies (/o = /n 

and therefore ||w —< C€{\\g\\i 2 (^Q'^ + ||s'||Lip)- Note that the computational 

complexity of subsequent solves to obtain is 0(e“'^(ln i)'^'*'^). 


Total computational complexity if a is periodic or ergodic with mixing length 

HP and e « with k > p. Under the assumptions of classical homogenization 
or HMM [33] (e.g. a is of period HP or a is ergodic with HP as mixing length), if 
the sets are chosen to match the period of a and the domain is rescaled so that 
1/H is an integer, then the entries of are invariant under periodic translations (or 
stationary if the medium is ergodic). Therefore, under these assumptions, as in classical 
homogenization, it sufficient to limit the computation of gamblets to periodicity cells (or 


51 



ergodicity cells with a tight control on mixing as in [47]). The resulting cost of obtaining 
^(fc),hom ^ W ^^2 ^ C C Q W 1^2 T 11 ^ 11 L i p ) ) is 

C>(iVln3'='A^i7P + e-‘^)ln'^+^ i). 
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